<<<<<<< HEAD COVID19 GER

Prepare data

Read and format prevalence data

df_ger_prev <- read_csv('GER_prevalence.csv')

df_ger_prev <- df_ger_prev %>%
  mutate(date = as.Date(date, "%d%b%Y"),
         kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day) %>%
  filter(date < '2020-04-01')

df_ger_prev %>% write_csv('ger_prev_kreis.csv')
df_ger_death <- read_csv('GER_prevalence.csv')

df_ger_death <-  df_ger_death %>%
  mutate(date = as.Date(date, "%d%b%Y"), 
         kreis = as.character(kreis))

df_ger_death_0930 <- df_ger_death %>% 
  filter(date == '2020-09-30') %>% 
  rename(case_rate_0930 = rate_day,
         death_rate_0930 = death_day) %>% 
  select(kreis, case_rate_0930, death_rate_0930)

df_ger_death_0331 <- df_ger_death %>% 
  filter(date == '2020-03-31') %>% 
  rename(case_rate_0331 = rate_day,
         death_rate_0331 = death_day) %>% 
  select(kreis, case_rate_0331, death_rate_0331)


df_ger_death <- merge(df_ger_death_0331, df_ger_death_0930)

# save df
df_ger_death %>% write_csv('ger_death_kreis_maps.csv')

Read and format scoial distancing data

fb_files <- list.files('../FB Data/GER individual files',
                       '*.csv', full.names = T)

df_ger_socdist <- fb_files %>% 
  map(read.csv) %>% bind_rows()

kreis_match <- read_csv('kreismatch.csv')

df_ger_socdist <- df_ger_socdist %>%
  select(-polygon_name) %>%
  plyr::join(kreis_match, by = 'polygon_id') %>%
  select(kreis, ds, all_day_bing_tiles_visited_relative_change,
         all_day_ratio_single_tile_users) %>%
  rename(date = ds,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(kreis = as.character(kreis),
         date = as.Date(date)) %>%
  arrange(kreis, date)

Read and format personality data

df_ger_pers <- read_csv('GER_personality_2.csv')

df_ger_pers %>% filter(frequ >= 50) %>% .$frequ %>% mean()
## [1] 284.0416
df_ger_pers <- df_ger_pers %>% 
  filter(frequ >= 50) %>%
  select(kreis, open, sci, extra, agree, neuro) %>%
  dplyr::rename(pers_o = open,
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct() %>%
  mutate(kreis = as.character(kreis))

Read and format county level controls

df_ger_ctrl <- read_csv('GER_controls_new.csv')

df_ger_ctrl <- df_ger_ctrl %>%
  select(-gdp, -cdu, -kreis_nme) %>%
  rename(tourism = tourism_beds,
         healthcare = hospital_beds,
         airport_dist = airport,
         conservative = afd, 
         popdens = popdens_new) %>%
  mutate(kreis = as.character(kreis))


df_ger_temp <- read_csv("GER_controls_weather_2.csv") %>% 
  select(kreis, kr_tamm_202003) %>% 
  rename(temp = kr_tamm_202003) %>%
  mutate(kreis = as.character(as.numeric(kreis)))

df_ger_ctrl <- df_ger_ctrl %>% left_join(df_ger_temp, 'kreis')

df_ger_pers %>% merge(df_ger_ctrl, by='kreis') %>% 
  write_csv('df_ger_pers_kreis.csv')

Merge prevalence data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge prevalence data
df_ger_prev <- df_ger_prev %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>%
  merge(df_dates, by='date') %>%
  arrange(kreis, date) %>%
  drop_na()

df_ger_prev %>% select(kreis) %>% distinct() %>% nrow()
## [1] 383
# join data frames
df_ger_death <- df_ger_death %>%
  plyr::join(df_ger_ctrl, by='kreis') %>%
  plyr::join(df_ger_pers, by='kreis') %>%
  drop_na()

# save df
df_ger_death %>%
  #select(kreis, case_rate, death_rate) %>%
  write_csv('ger_death_kreis_controls.csv')

Merge social distancing data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')

# merge socdist data
df_ger_socdist <- df_ger_socdist %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>% 
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_socdist %>% select(kreis) %>% distinct() %>% nrow()
## [1] 379

Control for weekend effect

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)

df_ger_loess <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_loess_2 <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_tiles ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  rename(loess_2 = loess) %>%
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_socdist <- df_ger_socdist %>% 
  merge(df_ger_loess, by=c('kreis', 'time')) %>% 
  merge(df_ger_loess_2, by=c('kreis', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile),
         socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess_2, socdist_tiles)) %>%
  arrange(kreis, time) %>% 
  select(-weekday)

df_ger_socdist <- df_ger_socdist %>% drop_na() %>% mutate(time = time-1)

Include cleaned social distancing data

df_ger_socdist <- df_ger_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

Define outcomes of interest

COVID - Extract onset

# get onset day
df_ger_onset_prev <- df_ger_prev %>% 
  group_by(kreis) %>% 
  filter(rate_day > 0.1) %>%
  summarise(onset_prev = min(time))
  
# merge with county data
df_ger_onset_prev <- df_ger_prev %>% 
  select(-date, -time, -rate_day) %>%
  distinct() %>% 
  left_join(df_ger_onset_prev, by = 'kreis')

# handle censored data
df_ger_onset_prev <- df_ger_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_ger_prev$date)))+1))

COVID - Extract slopes

# cut time series
df_ger_prev <- df_ger_prev %>% 
  filter(date > '2020-03-01' & date < '2020-04-01')

# extract slope prevalence
df_ger_slope_prev <- df_ger_prev %>% 
  split(.$kreis) %>% 
  map(~ lm(log(rate_day+1) ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>%
  rename(slope_prev = '.')

# merge with county data
df_ger_slope_prev <- df_ger_onset_prev %>% 
  left_join(df_ger_slope_prev, by = 'kreis')

# get unscaled object for descriptives
df_ger_prev_desc <- df_ger_slope_prev
# cut time series before onset
df_ger_slope_var <- df_ger_prev %>% 
  filter(date <= '2020-04-30') %>%
  group_by(kreis) %>% 
  filter(rate_day > 0.1) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30)

# extract slope prevalence
df_ger_slope_var <- df_ger_slope_var %>% 
  split(.$kreis) %>% 
  map(~ lm(log(rate_day+1) ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>% 
  rename(slope_prev_var = '.')

# merge with county data
df_ger_slope_prev <- df_ger_slope_prev %>% 
  left_join(df_ger_slope_var, by = 'kreis') %>%
  mutate(slope_prev_var = replace_na(slope_prev_var, 0))

Social Distancing - Change point analysis

# keep only counties with full data
kreis_complete <- df_ger_socdist %>% 
  group_by(kreis) %>% 
  summarise(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$kreis

# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist %>% 
  select(kreis, socdist_single_tile) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

df_ger_socdist_cpt_results_2 <- df_ger_socdist %>% 
  select(kreis, socdist_tiles) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_tiles),
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('kreis')

df_ger_socdist_cpt_day_2 <- df_ger_socdist_cpt_results_2 %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist_2 = '.') %>%
  rownames_to_column('kreis')

# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('kreis')

df_ger_socdist_cpt_mean_diff_2 <- df_ger_socdist_cpt_results_2 %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ -.[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist_2 = '.') %>%
  rownames_to_column('kreis')

# calculate means 
df_ger_socdist_mean <- df_ger_socdist %>%
  group_by(kreis) %>%
  summarise(mean_socdist = mean(socdist_single_tile))

df_ger_socdist_mean_2 <- df_ger_socdist %>%
  group_by(kreis) %>%
  summarise(mean_socdist_2 = -mean(socdist_tiles))

# merge with county data
df_ger_cpt_socdist <- df_ger_socdist %>%
  select(-date, -time, -socdist_single_tile, -socdist_tiles) %>%
  distinct() %>%
  left_join(df_ger_socdist_cpt_day, by='kreis') %>%
  left_join(df_ger_socdist_cpt_day_2, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff_2, by='kreis') %>%
  left_join(df_ger_socdist_mean, by='kreis') %>%
  left_join(df_ger_socdist_mean_2, by='kreis') %>%
  left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
  left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis')

# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), 
                                  as.numeric(diff(range(df_ger_socdist$date))), 
                                  cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 
                          as.numeric(diff(range(df_ger_socdist$date))), 0, 1))

Remove incomplete cases

df_ger_prev_unscaled <- inner_join(df_ger_slope_prev, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_death_unscaled <- inner_join(df_ger_death, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_socdist_unscaled <- inner_join(df_ger_cpt_socdist, df_ger_slope_prev %>% select(kreis), by='kreis')

kreis_list <- df_ger_prev_unscaled %>% select(kreis)

kreis_list %>% write_csv('kreis_list.csv')

df_ger_prev_unscaled %>% merge(df_ger_death_unscaled %>% select(kreis:death_rate_0930))  %>% 
  merge(df_ger_socdist_unscaled %>% select(kreis, cpt_day_socdist_2, mean_diff_socdist_2)) %>% 
  select(-event) %>% write_csv('ger_all_variables.csv')

Rescale Data

df_ger_prev_scaled <- df_ger_prev_unscaled %>% 
  mutate_at(vars(-kreis, -event), scale)
df_ger_death_scaled <- df_ger_death_unscaled %>% 
  mutate_at(vars(-kreis), scale)
df_ger_socdist_scaled <- df_ger_socdist_unscaled %>%
  mutate_at(vars(-kreis, -event), scale)

Calculate spatial lags

# create UDF to calculate lagged variables
calc_lags <- function(df, weights, cols){
  
  cols_only <- df %>% select(cols)
  cols_lag <- weights %*% as.matrix(cols_only) %>% 
  as.matrix() %>% as.data.frame()
names(cols_lag) <- paste0(names(cols_lag), '_lag')

return(cols_lag)
}
# read_weight matrix 
weight_mat_norm <- read_csv('GER_spatial_weights.csv')

weight_mat_norm <- weight_mat_norm %>% select(-kreis) %>% as.matrix()

dim(weight_mat_norm)
## [1] 379 379
# generate spatially lagged y 
y_only <- df_ger_prev_scaled %>% select(onset_prev,slope_prev,slope_prev_var) %>% names()
y_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_prev_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_prev_scaled <- cbind(df_ger_prev_scaled, y_lag, X_lag)
# generate spatially lagged y 
y_only <- df_ger_death_scaled %>% select(case_rate_0331:death_rate_0930) %>% names()
y_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_death_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_death_scaled <- cbind(df_ger_death_scaled, y_lag, X_lag)
# generate spatially lagged y 
y_only <- df_ger_socdist_scaled %>% select(cpt_day_socdist_2,mean_diff_socdist_2) %>% names()
y_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_socdist_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_socdist_scaled <- cbind(df_ger_socdist_scaled, y_lag, X_lag)
write_csv(df_ger_prev_scaled, 'df_ger_slope_prev.csv')
write_csv(df_ger_death_scaled, 'df_ger_death_scaled.csv')
write_csv(df_ger_socdist_scaled, 'df_ger_cpt_socdist.csv')

Run Specification Curve Analysis

# define function to calculate specification curve analyses for all traits
spec_calculate <- function(df, y, model, controls, all.comb = T){
  
  spec_results_o <- run_specs(df = df,
                       y = c(y), x = c("pers_o"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_c <- run_specs(df = df, 
                       y = c(y), x = c("pers_c"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_e <- run_specs(df = df, 
                       y = c(y), x = c("pers_e"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_a <- run_specs(df = df, 
                       y = c(y), x = c("pers_a"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_n <- run_specs(df = df, 
                       y = c(y), x = c("pers_n"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results <- list(spec_results_o, spec_results_c, spec_results_e, 
                      spec_results_a, spec_results_n)
  
  names(spec_results) <- list('spec_results_o', 'spec_results_c', 'spec_results_e', 
                      'spec_results_a', 'spec_results_n')

  return(spec_results)
  }
# adapted from specr package code
format_results <- function(df, var, null = 0, desc = FALSE) {

  # rank specs
  if (isFALSE(desc)) {
    df <- df %>%
      dplyr::arrange(!! var)
  } else {
    df <- df %>%
      dplyr::arrange(desc(!! var))
  }

  # create rank variable and color significance
  df <- df %>%
    dplyr::mutate(specifications = 1:nrow(df),
                  color = case_when(conf.low > null ~ "#377eb8",
                                    conf.high < null ~ "#e41a1c",
                                    TRUE ~ "darkgrey"))
  return(df)
}


# define function to plot single specification curve
plot_curves <- function(results, filename, hr=F){
  
  file_path_eps <- paste0("../../Plots/GER/", filename,".eps")
  file_path_pdf <- paste0("../../Plots/GER/", filename,".pdf")
  file_path_png <- paste0("../../Plots/GER/", filename,".png")

  if(hr==F){
    
    results %>%
    format_results(var = .$estimate, null = 0, desc = F) %>%
    ggplot(aes(x = specifications,
               y = estimate,
               ymin = conf.low,
               ymax = conf.high,
               color = color)) +
    geom_point(aes(color = color),
               size = 1) +
    theme_minimal() +
    scale_color_identity() +
    theme(strip.text = element_blank(),
          axis.line = element_line("black", size = .5),
          legend.position = "none",
          panel.spacing = unit(.75, "lines"),
          axis.text = element_text(colour = "black")) +
    labs(x = "") +
      geom_pointrange(alpha = 0.05,
                      size = .6,
                      fatten = 1) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    labs(x = "", y = "standarized coefficient") + 
      coord_fixed(ratio = 2000, ylim = c(-0.4, 0.4)) +
      theme(panel.border = element_blank(), panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            axis.line = element_line(colour = "black")) +
      ggsave(file=file_path_eps, device = 'eps')+
      ggsave(file=file_path_pdf, device = 'pdf')+
      ggsave(file=file_path_png, device = 'png')
    
  }else{
    
    results %>%
    format_results(var = .$estimate, null = 1, desc = F) %>%
    ggplot(aes(x = specifications,
               y = estimate,
               ymin = conf.low,
               ymax = conf.high,
               color = color)) +
    geom_point(aes(color = color),
               size = 1) +
    theme_minimal() +
    scale_color_identity() +
    theme(strip.text = element_blank(),
          axis.line = element_line("black", size = .5),
          legend.position = "none",
          panel.spacing = unit(.75, "lines"),
          axis.text = element_text(colour = "black")) +
    labs(x = "") +
      geom_pointrange(alpha = 0.05,
                      size = .6,
                      fatten = 1) +
    geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
    labs(x = "", y = "standarized coefficient") + 
      coord_fixed(ratio = 2000, ylim = c(0.6, 1.4)) +
      theme(panel.border = element_blank(), panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(), 
      axis.line = element_line(colour = "black")) +
      ggsave(file=file_path_eps, device = 'eps')+
      ggsave(file=file_path_pdf, device = 'pdf')+
      ggsave(file=file_path_png, device = 'png')  
    }
  
}

# define function to plot all curves in list
plot_all_curves <- function(ls, analysis, hr=F){
  
  pers <- c('o', 'c', 'e', 'a', 'n')
  filenames <- as.list(paste0(analysis, '_', pers))
  hr <- as.list(rep(hr, 5))
  pmap(list(ls, filenames, hr), plot_curves)
}
# define function to calculate summary stats for single specification curve
calc_summary <- function(df){
  
  dft <- df %>% select(estimate:fit_nobs)
  dft <- dft %>% mutate(significant = as.numeric(p.value<0.05),
                        positive = as.numeric(statistic>0), 
                        negative = as.numeric(statistic<0), 
                        significant_positive = as.numeric(p.value<0.05 & statistic>0),
                        significant_negative = as.numeric(p.value<0.05 & statistic<0))
  
  mean_temp <- dft %>% map_if(is.numeric, mean, .else=NULL) %>% as.data.frame()
  sd_temp <- dft %>% map_if(is.numeric, sd, .else=NULL) %>% as.data.frame()
  
  df_temp <- rbind(mean_temp, sd_temp)
  row.names(df_temp) <- c('mean', 'sd')
  
  return(df_temp)
}

# define function to calculate summary stats for all curves in list
calc_all_summaries <- function(ls){
  
  ls <- ls %>% map(calc_summary)
  names(ls) <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
  return(ls)
}
coef_to_hr <- function(df){
  
  df %>% mutate(conf.low = exp(estimate-1.96*std.error),
           conf.high = exp(estimate+1.96*std.error),
           estimate = exp(estimate)) 
}

filter_socdist <- function(df){
  
  df %>% filter(str_detect(controls, 'onset_prev') & 
                  str_detect(controls, 'slope_prev'))
}
covariates <- c("airport_dist", "conservative", 'male', 'age', 'popdens',
                'manufact', 'tourism', 'academics', 'medinc', 'healthcare',
                'temp')

Predict Onset

cox_model <- function(formula, data){
  formula <- as.formula(formula)
  coxph(formula = formula, data = data)}
cox_onset_prev <- spec_calculate(df = df_ger_prev_scaled, 
               y = "Surv(onset_prev, event)", 
               model = "cox_model", 
               controls = covariates %>% append('onset_prev_lag'),
               all.comb = T)

cox_onset_prev_hr <- cox_onset_prev %>% map(coef_to_hr)

calc_all_summaries(cox_onset_prev_hr)
## $pers_o
##        estimate   std.error statistic   p.value   conf.low  conf.high fit_n
## mean 1.02726691 0.055320547 0.4779412 0.5231112 0.92180104 1.14483537   379
## sd   0.05390301 0.002840416 0.9535330 0.2800335 0.05018495 0.05845821     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         131.39127    3.811112e-06         140.8145
## sd            0          36.26568    5.698238e-05          41.9226
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   4.156490e-06          137.52997     4.234893e-06                   NA
## sd     5.770949e-05           40.56007     5.819916e-05                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.2895930         0.9999496       0.7195157
## sd                   NA     0.0722395         0.0000000       0.0381976
##      fit_std.error.concordance  fit_logLik   fit_AIC    fit_BIC fit_nobs
## mean               0.015635755 -1809.51851 3633.0370 3660.59978      379
## sd                 0.001317056    18.13284   34.3554   31.38233        0
##      significant  positive  negative significant_positive significant_negative
## mean  0.08911133 0.6496582 0.3503418           0.08911133                    0
## sd    0.28493915 0.4771352 0.4771352           0.28493915                    0
## 
## $pers_c
##        estimate   std.error  statistic   p.value   conf.low  conf.high fit_n
## mean 0.94922252 0.053595445 -0.9954995 0.3479716 0.85456571 1.05437271   379
## sd   0.04503002 0.001467328  0.8933022 0.2887067 0.04054931 0.05018395     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         132.01667    0.0001954177         140.8434
## sd            0          37.85635    0.0029957287          43.1195
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   0.0001852582          137.59893     0.0001881611                   NA
## sd     0.0028016025           41.78484     0.0028121428                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.29044877         0.9999496      0.71826746
## sd                   NA    0.07571748         0.0000000      0.04121107
##      fit_std.error.concordance  fit_logLik    fit_AIC   fit_BIC fit_nobs
## mean               0.015744358 -1809.20581 3632.41162 3659.9744      379
## sd                 0.001364156    18.92817   35.97484   33.0281        0
##      significant  positive  negative significant_positive significant_negative
## mean   0.1445312 0.1584473 0.8415527                    0            0.1445312
## sd     0.3516705 0.3652045 0.3652045                    0            0.3516705
## 
## $pers_e
##        estimate   std.error statistic   p.value   conf.low  conf.high fit_n
## mean 1.10484855 0.052014824 1.9131109 0.1418816 0.99784541 1.22333370   379
## sd   0.05117709 0.001260917 0.9230424 0.1475011 0.04805772 0.05449954     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         134.80506    6.667657e-08        144.44569
## sd            0          33.69771    7.667802e-07         39.58589
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   1.086854e-07          141.00631     1.184551e-07                   NA
## sd     1.098564e-06           38.16651     1.180765e-06                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.29642439         0.9999496      0.72680567
## sd                   NA    0.06635634         0.0000000      0.03214107
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.015317233 -1807.81162 3629.62324 3657.18599      379
## sd                 0.001142571    16.84885   31.78596   28.87968        0
##      significant positive negative significant_positive significant_negative
## mean   0.3879395        1        0            0.3879395                    0
## sd     0.4873401        0        0            0.4873401                    0
## 
## $pers_a
##      estimate   std.error statistic  p.value   conf.low  conf.high fit_n
## mean 1.006651 0.053837899 0.1292402 0.645567 0.90592390 1.11859211   379
## sd   0.032338 0.001749636 0.6071196 0.238925 0.03165915 0.03294262     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          130.6175    0.0006586847        140.37322
## sd            0           37.6386    0.0123791555         43.20648
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   0.0006630851          137.15705      0.000668293                   NA
## sd     0.0122500086           41.79493      0.012237339                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.2878630         0.9999496      0.71879483
## sd                   NA     0.0756924         0.0000000      0.04050933
##      fit_std.error.concordance fit_logLik   fit_AIC    fit_BIC fit_nobs
## mean               0.015616239 -1809.9054 3633.8108 3661.37351      379
## sd                 0.001276486    18.8193   35.7401   32.76038        0
##      significant  positive  negative significant_positive significant_negative
## mean           0 0.4870605 0.5129395                    0                    0
## sd             0 0.4998936 0.4998936                    0                    0
## 
## $pers_n
##       estimate    std.error  statistic    p.value   conf.low  conf.high fit_n
## mean 0.8618104 0.0536108475 -2.7813176 0.01334193 0.77585634 0.95728861   379
## sd   0.0256510 0.0006767025  0.5515346 0.01683804 0.02333785 0.02824371     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         138.42138    8.297546e-07        148.34185
## sd            0          35.84265    1.236878e-05         41.64845
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   9.623814e-07          144.57863     1.004892e-06                   NA
## sd     1.340854e-05           40.36034     1.343165e-05                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.30271384         0.9999496      0.72815172
## sd                   NA    0.07048457         0.0000000      0.03346175
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.015109031 -1806.00346 3626.00692 3653.56967      379
## sd                 0.001124059    17.92132   33.95795   31.05319        0
##      significant positive negative significant_positive significant_negative
## mean   0.9492188        0        1                    0            0.9492188
## sd     0.2195775        0        0                    0            0.2195775
plot_all_curves(cox_onset_prev_hr, 'cox_onset_prev_hr', hr = T)
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev_lag, 
                             data = df_ger_prev_scaled)

summary(cox_onset_prev_ctrl)
## Call:
## coxph(formula = Surv(onset_prev, event) ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + onset_prev_lag, data = df_ger_prev_scaled)
## 
##   n= 379, number of events= 379 
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)    
## airport_dist   -0.09147   0.91259  0.07045 -1.298 0.194133    
## conservative   -0.08051   0.92264  0.08006 -1.006 0.314608    
## male           -0.13923   0.87003  0.07451 -1.869 0.061664 .  
## age            -0.25167   0.77750  0.09118 -2.760 0.005779 ** 
## popdens        -0.09156   0.91251  0.07416 -1.235 0.216965    
## manufact        0.00861   1.00865  0.05321  0.162 0.871458    
## tourism        -0.11689   0.88968  0.05988 -1.952 0.050933 .  
## academics       0.07279   1.07551  0.08113  0.897 0.369583    
## medinc          0.29682   1.34558  0.07898  3.758 0.000171 ***
## healthcare     -0.06525   0.93683  0.06938 -0.940 0.347005    
## temp           -0.07546   0.92731  0.07139 -1.057 0.290488    
## onset_prev_lag -0.36945   0.69111  0.08959 -4.124 3.73e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## airport_dist      0.9126     1.0958    0.7949    1.0477
## conservative      0.9226     1.0838    0.7886    1.0794
## male              0.8700     1.1494    0.7518    1.0068
## age               0.7775     1.2862    0.6503    0.9296
## popdens           0.9125     1.0959    0.7891    1.0553
## manufact          1.0086     0.9914    0.9087    1.1195
## tourism           0.8897     1.1240    0.7912    1.0005
## academics         1.0755     0.9298    0.9174    1.2609
## medinc            1.3456     0.7432    1.1526    1.5709
## healthcare        0.9368     1.0674    0.8177    1.0733
## temp              0.9273     1.0784    0.8062    1.0666
## onset_prev_lag    0.6911     1.4469    0.5798    0.8238
## 
## Concordance= 0.755  (se = 0.015 )
## Likelihood ratio test= 177.8  on 12 df,   p=<2e-16
## Wald test            = 183.9  on 12 df,   p=<2e-16
## Score (logrank) test = 193.2  on 12 df,   p=<2e-16

Predict Slopes

# fixed time windows
lm_slope_prev <- spec_calculate(df = df_ger_prev_scaled, 
               y = "slope_prev", 
               model = "lm", 
               controls = covariates %>% append('slope_prev_lag'),
               all.comb = T)

calc_all_summaries(lm_slope_prev)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04770701 0.047105679 0.9732636 0.4496175 -0.04492057 0.14033459
## sd   0.04866155 0.003961232 0.9445598 0.3249956  0.04463451 0.05352556
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3762331         0.3647529 0.79419554      34.39651 0.0000103227
## sd       0.1110695         0.1112149 0.06709402      12.97273 0.0001987136
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.11260 908.22521 943.66303    235.78389      371.000000
## sd   1.732262   31.40818  60.85522  57.41092     41.98425        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1411133 0.8693848 0.1306152            0.1411133
## sd          0   0.3481809 0.3370202 0.3370202            0.3481809
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.02472864 0.042382347 -0.5866914 0.4324615 -0.10806835 0.05861107
## sd    0.04143443 0.003375866  0.9654986 0.2623563  0.04195917 0.04196624
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean     0.3747063         0.3632077 0.79495961      34.26618 0.001065338
## sd       0.1157000         0.1159113 0.06951753      13.36336 0.018659191
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.39270 908.78539 944.22322    236.36101      371.000000
## sd   1.732262   32.35788  62.76145  59.31843     43.73461        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.09521484 0.3103027 0.6896973                    0
## sd          0  0.29354730 0.4626740 0.4626740                    0
##      significant_negative
## mean           0.09521484
## sd             0.29354730
## 
## $pers_e
##        estimate   std.error statistic   p.value   conf.low  conf.high
## mean 0.09475842 0.042916475  2.145207 0.1152644 0.01036841 0.17914843
## sd   0.05415151 0.002851794  1.072831 0.1384433 0.04912212 0.05928462
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3835310         0.3721721 0.78983225      35.33176 6.972021e-08
## sd       0.1038293         0.1038635 0.06319739      12.45035 1.032508e-06
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -443.15352 904.30703 939.74486    233.02529      371.000000
## sd   1.732262   29.85157  57.72422  54.26548     39.24747        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.4655762        1        0            0.4655762
## sd          0   0.4988745        0        0            0.4988745
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04444294 0.042531589 1.0138143 0.4314829 -0.03919023 0.12807612
## sd   0.03737716 0.003433074 0.8010556 0.3214330  0.03396617 0.04161167
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean     0.3757420         0.3642556 0.7944327      34.36378 0.0001083107
## sd       0.1128081         0.1129825 0.0679861      13.13045 0.0019291217
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.19538 908.3908 943.82859    235.96953      371.000000
## sd   1.732262   31.75005  61.5464  58.11313     42.64145        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379   0.1437988 0.9528809 0.04711914            0.1437988
## sd          0   0.3509285 0.2119195 0.21191954            0.3509285
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error statistic   p.value    conf.low  conf.high
## mean -0.06378630 0.042452946 -1.458345 0.3027720 -0.14726483 0.01969223
## sd    0.04686406 0.003144902  0.979003 0.2759218  0.05132561 0.04283259
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3785076         0.3670643 0.79281621      34.71477 2.890617e-06
## sd       0.1093210         0.1094544 0.06617581      12.90910 5.052404e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -444.48386 906.96772 942.40555    234.92414      371.000000
## sd   1.732262   31.04716  60.14311  56.72546     41.32336        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.3442383        0        1                    0
## sd          0   0.4751772        0        0                    0
##      significant_negative
## mean            0.3442383
## sd              0.4751772
plot_all_curves(lm_slope_prev, 'lm_slope_prev')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_slope_prev_ctrl <- lm(slope_prev ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + slope_prev_lag, 
                             data = df_ger_prev_scaled)

summary(lm_slope_prev_ctrl)
## 
## Call:
## lm(formula = slope_prev ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + slope_prev_lag, data = df_ger_prev_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5022 -0.4555 -0.1231  0.3315  3.8222 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -0.007574   0.036347  -0.208 0.835056    
## airport_dist   -0.010112   0.046450  -0.218 0.827784    
## conservative    0.007307   0.056396   0.130 0.896981    
## male           -0.085540   0.053273  -1.606 0.109207    
## age            -0.296574   0.069586  -4.262 2.58e-05 ***
## popdens        -0.116212   0.058141  -1.999 0.046368 *  
## manufact        0.061335   0.042127   1.456 0.146268    
## tourism        -0.126881   0.041916  -3.027 0.002645 ** 
## academics       0.025668   0.056486   0.454 0.649799    
## medinc          0.239428   0.052239   4.583 6.29e-06 ***
## healthcare     -0.018543   0.047757  -0.388 0.698034    
## temp           -0.174761   0.049146  -3.556 0.000426 ***
## slope_prev_lag  0.432857   0.060688   7.133 5.31e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7073 on 366 degrees of freedom
## Multiple R-squared:  0.5156, Adjusted R-squared:  0.4997 
## F-statistic: 32.47 on 12 and 366 DF,  p-value: < 2.2e-16
# variable time windows
lm_slope_prev_var <- spec_calculate(df = df_ger_prev_scaled, 
               y = "slope_prev_var", 
               model = "lm", 
               controls = covariates %>% append('slope_prev_var_lag'),
               all.comb = T)

calc_all_summaries(lm_slope_prev_var)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low conf.high
## mean 0.03138877 0.048027782 0.6168286 0.5240564 -0.06305202 0.1258296
## sd   0.04551793 0.003904359 0.8881225 0.3037374  0.04165726 0.0502624
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3513107         0.3393623 0.81020283      30.71510 0.0001379822
## sd       0.1089591         0.1089781 0.06488532      11.75029 0.0027232807
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.79671 923.59341 959.03124    245.20455      371.000000
## sd   1.732262   29.99134  57.97282  54.44501     41.18655        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.07739258 0.7424316 0.2575684           0.07739258
## sd          0  0.26724596 0.4373484 0.4373484           0.26724596
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.005412083 0.043251658 -0.1122928 0.4938687 -0.09046118 0.07963702
## sd    0.041568173 0.003203958  0.9510700 0.2765537  0.04319003 0.04086338
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean     0.3505259         0.3385690 0.81057266      30.65762 0.001475911
## sd       0.1118288         0.1118893 0.06636315      11.99597 0.025347073
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.92184 923.8437 959.28150    245.50122      371.000000
## sd   1.732262   30.55817  59.1121  55.58766     42.27129        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.04272461 0.4780273 0.5219727                    0
## sd          0  0.20226024 0.4995780 0.4995780                    0
##      significant_negative
## mean           0.04272461
## sd             0.20226024
## 
## $pers_e
##        estimate   std.error statistic   p.value    conf.low conf.high
## mean 0.06672988 0.043867749  1.464399 0.3014781 -0.01953070 0.1529904
## sd   0.05252357 0.002821869  1.068661 0.2892869  0.04787854 0.0573294
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3555855         0.3437053 0.80774108      31.20591 2.225108e-06
## sd       0.1036970         0.1036339 0.06204795      11.33341 3.180183e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -451.73375 921.46749 956.90532    243.58870      371.000000
## sd   1.732262   28.85099  55.67972  52.14376     39.19745        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379   0.2744141 0.9606934 0.03930664            0.2744141
## sd          0   0.4462730 0.1943472 0.19434724            0.4462730
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05484896 0.043310646 1.2488281 0.3161042 -0.03031614 0.14001406
## sd   0.03452171 0.003320092 0.7388706 0.2662508  0.03219300 0.03784622
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean     0.3528556         0.3409354 0.8092101      30.95031 0.0001543231
## sd       0.1093621         0.1094004 0.0651510      11.88643 0.0025729140
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000  -452.3210 922.64193 958.07975    244.62057      371.000000
## sd   1.732262    30.1185  58.24095  54.74094     41.33889        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1984863 0.9831543 0.0168457            0.1984863
## sd          0   0.3989090 0.1287089 0.1287089            0.3989090
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.01945947 0.043371756 -0.4091114 0.4680019 -0.10474473 0.06582579
## sd    0.04490624 0.003125042  0.9702821 0.2636699  0.04887248 0.04147441
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3511661         0.3392152 0.81027906      30.70800 7.639954e-05
## sd       0.1092723         0.1093028 0.06506678      11.78065 1.248112e-03
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000  -452.8261 923.65224 959.09007    245.25920      371.000000
## sd   1.732262    30.0687  58.13108  54.60971     41.30493        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.08203125 0.4973145 0.5026855                    0
## sd          0  0.27444583 0.5000538 0.5000538                    0
##      significant_negative
## mean           0.08203125
## sd             0.27444583
plot_all_curves(lm_slope_prev_var, 'lm_slope_prev_var')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_slope_prev_var_ctrl <- lm(slope_prev_var ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + slope_prev_var_lag, 
                             data = df_ger_prev_scaled)

summary(lm_slope_prev_var_ctrl)
## 
## Call:
## lm(formula = slope_prev_var ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + slope_prev_var_lag, data = df_ger_prev_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4605 -0.4684 -0.0894  0.3536  3.8629 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.0107049  0.0370087  -0.289 0.772552    
## airport_dist       -0.0006466  0.0472676  -0.014 0.989094    
## conservative        0.0578272  0.0574076   1.007 0.314452    
## male               -0.1063530  0.0542440  -1.961 0.050679 .  
## age                -0.3348173  0.0704594  -4.752 2.90e-06 ***
## popdens            -0.1144899  0.0591642  -1.935 0.053745 .  
## manufact            0.0844884  0.0428457   1.972 0.049371 *  
## tourism            -0.1367335  0.0426814  -3.204 0.001476 ** 
## academics          -0.0251356  0.0576842  -0.436 0.663279    
## medinc              0.2276969  0.0521603   4.365 1.65e-05 ***
## healthcare         -0.0223008  0.0487212  -0.458 0.647423    
## temp               -0.1940683  0.0505286  -3.841 0.000145 ***
## slope_prev_var_lag  0.4329906  0.0602916   7.182 3.88e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7199 on 366 degrees of freedom
## Multiple R-squared:  0.4982, Adjusted R-squared:  0.4817 
## F-statistic: 30.28 on 12 and 366 DF,  p-value: < 2.2e-16

Predict Cases

lm_cases_0331 <- spec_calculate(df = df_ger_death_scaled, 
               y = "case_rate_0331", 
               model = "lm", 
               controls = covariates %>% append('case_rate_0331_lag'),
               all.comb = T)

calc_all_summaries(lm_cases_0331)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.03158699 0.051022780 0.5994671 0.5627758 -0.06874312 0.13191710
## sd   0.04075752 0.003039436 0.7655479 0.2958935  0.03824516 0.04394422
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27125522         0.2577158 0.86029794     20.714339 0.0002426819
## sd      0.08294606         0.0823706 0.04660655      7.515557 0.0038974652
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.16370 970.32739 1005.76522    275.46553      371.000000
## sd   1.732262   20.54957  39.18622   36.10197     31.35361        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.06103516 0.7888184 0.2111816           0.06103516
## sd          0  0.23942402 0.4081968 0.4081968           0.23942402
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.004091242 0.045901847 -0.08596589 0.5796038 -0.09435164 0.08616916
## sd    0.034700872 0.002269187  0.74910889 0.2613804  0.03536478 0.03460409
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean    0.27028449        0.25673269 0.8607943     20.626725 0.002355802
## sd      0.08560184        0.08506666 0.0479699      7.724115 0.031524022
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.35290 970.70581 1006.14364     275.8325      371.000000
## sd   1.732262   21.07044  40.23365   37.14286      32.3575        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.4645996 0.5354004                    0
## sd          0           0 0.4988061 0.4988061                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate   std.error statistic   p.value     conf.low  conf.high
## mean 0.08714459 0.046581747 1.8422270 0.1520497 -0.004452756 0.17874193
## sd   0.04647008 0.001763032 0.9112771 0.1493776  0.043558178 0.04945351
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27812366        0.26469914 0.85640497     21.418135 2.574047e-06
## sd      0.07660416        0.07592723 0.04326494      7.115188 3.140562e-05
##        fit_df fit_logLik  fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -474.51363 967.0273 1002.46509    272.86926      371.000000
## sd   1.732262   19.24243  36.5560   33.48926     28.95637        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.3618164        1        0            0.3618164
## sd          0   0.4805847        0        0            0.4805847
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.07275316 0.045946320 1.5732114 0.1733925 -0.01759469 0.16310101
## sd   0.02970329 0.002372876 0.5983738 0.1394271  0.02789450 0.03209363
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27494180        0.26146891 0.85812315     21.120691 0.0001372159
## sd      0.08246331        0.08190365 0.04643572      7.589597 0.0017659528
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.20596 968.41191 1003.8497    274.07200      371.000000
## sd   1.732262   20.51326  39.13057   36.0854     31.17113        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2634277        1        0            0.2634277
## sd          0   0.4405462        0        0            0.4405462
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error statistic   p.value    conf.low  conf.high
## mean -0.06299523 0.045987096 -1.352350 0.2837433 -0.15342327 0.02743280
## sd    0.03825584 0.002107881  0.770264 0.2261862  0.04073718 0.03608117
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27430368        0.26081464 0.85855349      21.04466 2.747966e-05
## sd      0.08078661        0.08020112 0.04551672       7.43662 3.526301e-04
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.41660 968.83321 1004.27103    274.31321      371.000000
## sd   1.732262   20.13405  38.37129   35.33896     30.53734        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2299805        0        1                    0
## sd          0   0.4208714        0        0                    0
##      significant_negative
## mean            0.2299805
## sd              0.4208714
plot_all_curves(lm_cases_0331, 'lm_cases_0331')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_cases_0331_ctrl <- lm(case_rate_0331 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + case_rate_0331_lag, 
                             data = df_ger_death_scaled)

summary(lm_cases_0331_ctrl)
## 
## Call:
## lm(formula = case_rate_0331 ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + case_rate_0331_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4908 -0.4246 -0.1322  0.1935  5.7140 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.005574   0.041192  -0.135 0.892434    
## airport_dist        0.007239   0.052639   0.138 0.890700    
## conservative        0.021561   0.064001   0.337 0.736398    
## male               -0.079280   0.060385  -1.313 0.190032    
## age                -0.287805   0.078147  -3.683 0.000265 ***
## popdens            -0.094136   0.065815  -1.430 0.153484    
## manufact            0.045203   0.047719   0.947 0.344117    
## tourism            -0.094058   0.047518  -1.979 0.048517 *  
## academics          -0.011510   0.063625  -0.181 0.856548    
## medinc              0.184173   0.059260   3.108 0.002032 ** 
## healthcare         -0.036733   0.054269  -0.677 0.498913    
## temp               -0.180424   0.055771  -3.235 0.001327 ** 
## case_rate_0331_lag  0.404243   0.071769   5.633 3.54e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8017 on 366 degrees of freedom
## Multiple R-squared:  0.3777, Adjusted R-squared:  0.3573 
## F-statistic: 18.51 on 12 and 366 DF,  p-value: < 2.2e-16
lm_cases_0930 <- spec_calculate(df = df_ger_death_scaled, 
               y = "case_rate_0930", 
               model = "lm", 
               controls = covariates %>% append('case_rate_0930_lag'),
               all.comb = T)

calc_all_summaries(lm_cases_0930)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04703511 0.048759556 0.9328083 0.4213025 -0.04884462 0.14291485
## sd   0.04623303 0.003708065 0.8870397 0.3049813  0.04263197 0.05063419
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3311692         0.3188131 0.82304322      28.08371 0.0000322385
## sd       0.1042954         0.1042904 0.06154394      11.21308 0.0004799432
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -458.89825 935.79650 971.23433    252.81805      371.000000
## sd   1.732262   28.15626  54.43081  51.21137     39.42368        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1430664 0.8413086 0.1586914            0.1430664
## sd          0   0.3501833 0.3654327 0.3654327            0.3501833
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.002610167 0.044018978 -0.02673773 0.5244260 -0.08916812 0.08394778
## sd    0.040295483 0.003040096  0.89565674 0.2664612  0.04373607 0.03749754
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3292839         0.3169004 0.82408293      27.87911 0.0006614882
## sd       0.1074214         0.1074516 0.06314984      11.42535 0.0139028188
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.32868 936.65736 972.09518    253.53069      371.000000
## sd   1.732262   28.77586  55.66624  52.42724     40.60527        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.0324707 0.5195312 0.4804688                    0
## sd          0   0.1772682 0.4996794 0.4996794                    0
##      significant_negative
## mean            0.0324707
## sd              0.1772682
## 
## $pers_e
##        estimate   std.error statistic    p.value   conf.low  conf.high
## mean 0.11873043 0.044310503 2.6373935 0.03352240 0.03159922 0.20586163
## sd   0.04761501 0.002634405 0.9060562 0.04271978 0.04339097 0.05201209
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.34275305        0.33059563 0.81616073      29.48154 1.913204e-08
## sd      0.09612161        0.09600844 0.05733096      10.80161 3.024722e-07
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -455.83276 929.6655 965.10334    248.43935      371.000000
## sd   1.732262   26.55416  51.2221  48.02528     36.33397        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.7509766        1        0            0.7509766
## sd          0   0.4325002        0        0            0.4325002
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean 0.001807263 0.044142706 0.02930601 0.6298433 -0.08499398 0.08860851
## sd   0.029309078 0.003273494 0.63816227 0.2379830  0.02826462 0.03165452
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean      0.328552         0.3161559 0.8245256      27.78624 0.0008432176
## sd        0.107648         0.1076773 0.0632661      11.42428 0.0168326235
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.52910 937.0582 972.49602    253.80733      371.000000
## sd   1.732262   28.81963  55.7506  52.50393     40.69093        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.5290527 0.4709473                    0
## sd          0           0 0.4992162 0.4992162                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.001840958 0.044147329 -0.01261312 0.5589962 -0.08865129 0.08496938
## sd    0.037168325 0.003113077  0.79081236 0.2465463  0.04075222 0.03430969
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3290563         0.3166673 0.82425385      27.83814 0.0001595196
## sd       0.1066349         0.1066540 0.06276337      11.33910 0.0028612546
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.41897 936.83794 972.27576    253.61671      371.000000
## sd   1.732262   28.63422  55.37859  52.13274     40.30798        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.01220703 0.5891113 0.4108887                    0
## sd          0  0.10982242 0.4920552 0.4920552                    0
##      significant_negative
## mean           0.01220703
## sd             0.10982242
plot_all_curves(lm_cases_0930, 'lm_cases_0930')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_cases_0930_ctrl <- lm(case_rate_0930 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + case_rate_0930_lag, 
                             data = df_ger_death_scaled)

summary(lm_cases_0930_ctrl)
## 
## Call:
## lm(formula = case_rate_0930 ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + case_rate_0930_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6305 -0.4216 -0.1215  0.3040  3.5774 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.003935   0.038008  -0.104 0.917604    
## airport_dist       -0.045205   0.048579  -0.931 0.352702    
## conservative        0.102075   0.058827   1.735 0.083551 .  
## male               -0.034565   0.055737  -0.620 0.535548    
## age                -0.350952   0.071996  -4.875 1.63e-06 ***
## popdens             0.152294   0.061336   2.483 0.013478 *  
## manufact            0.062070   0.044049   1.409 0.159651    
## tourism            -0.109885   0.043906  -2.503 0.012760 *  
## academics          -0.217018   0.059681  -3.636 0.000316 ***
## medinc              0.190831   0.056091   3.402 0.000742 ***
## healthcare          0.013822   0.049981   0.277 0.782285    
## temp               -0.156873   0.051285  -3.059 0.002385 ** 
## case_rate_0930_lag  0.479824   0.068459   7.009 1.16e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7399 on 366 degrees of freedom
## Multiple R-squared:   0.47,  Adjusted R-squared:  0.4526 
## F-statistic: 27.05 on 12 and 366 DF,  p-value: < 2.2e-16

Predict Deaths

lm_deaths_0331 <- spec_calculate(df = df_ger_death_scaled, 
               y = "death_rate_0331", 
               model = "lm", 
               controls = covariates %>% append('death_rate_0331_lag'),
               all.comb = T)

calc_all_summaries(lm_deaths_0331)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05926579 0.056344799  1.048407 0.3391156 -0.05152945 0.17006104
## sd   0.02649895 0.002223663  0.465681 0.2207861  0.02562865 0.02803226
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11230093        0.09564048 0.95084691      6.817236 0.003741489
## sd      0.03249807        0.03024746 0.01580279      1.831876 0.022251556
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.578479 1047.15696 1082.59478    335.55025      371.000000
## sd   1.732262    6.852099   11.52894    9.74724     12.28427        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379  0.02636719 0.9821777 0.01782227           0.02636719
## sd          0  0.16024428 0.1323212 0.13232123           0.16024428
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate    std.error  statistic   p.value    conf.low  conf.high
## mean -0.02176714 0.0507656910 -0.4272184 0.6489583 -0.12159172 0.07805745
## sd    0.02201730 0.0008507798  0.4327864 0.2340294  0.02246288 0.02169193
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11002111        0.09331891 0.95206101      6.651521  0.00963241
## sd      0.03320389        0.03098238 0.01615476      1.883209  0.05665206
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -515.059950 1048.11990 1083.557726    336.41202      371.000000
## sd   1.732262    6.973763   11.79201    9.998775     12.55107        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1816406 0.8183594                    0
## sd          0           0 0.3855951 0.3855951                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate    std.error statistic   p.value    conf.low  conf.high
## mean 0.06219816 0.0515241068 1.2102704 0.2768198 -0.03911776 0.16351408
## sd   0.02481655 0.0007406798 0.4899089 0.1810235  0.02544552 0.02425882
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11323239        0.09658630 0.95035178      6.898352 0.002436617
## sd      0.03213563        0.02994777 0.01566049      1.848994 0.014602120
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.381924 1046.76385 1082.201674    335.19815      371.000000
## sd   1.732262    6.787753   11.45944    9.857845     12.14727        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379  0.09594727        1        0           0.09594727
## sd          0  0.29455487        0        0           0.29455487
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.07596961 0.050837168 1.4956105 0.1615796 -0.02399553 0.17593475
## sd   0.01854741 0.001137139 0.3694743 0.1092908  0.01880931 0.01855325
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean    0.11477591        0.09816324 0.9495095      6.996946 0.003652699
## sd      0.03348276        0.03131897 0.0163866      1.927123 0.022385471
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.040955 1046.08191 1081.51974    334.61471      371.000000
## sd   1.732262    7.080364   12.03403   10.26139     12.65648        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.1005859        1        0            0.1005859
## sd          0   0.3008164        0        0            0.3008164
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.05389933 0.050903492 -1.0622345 0.3222672 -0.15399488 0.04619623
## sd    0.01914286 0.001035423  0.3829244 0.1659477  0.01845582 0.02001427
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11220076        0.09553655 0.95090074      6.820286 0.003463388
## sd      0.03252479        0.03033547 0.01585213      1.860572 0.020486709
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.599597 1047.19919 1082.637020    335.58811      371.000000
## sd   1.732262    6.860555   11.59847    9.943315     12.29437        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379 0.005859375        0        1                    0
## sd          0 0.076331286        0        0                    0
##      significant_negative
## mean          0.005859375
## sd            0.076331286
plot_all_curves(lm_deaths_0331, 'lm_deaths_0331')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_deaths_0331_ctrl <- lm(death_rate_0331 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + death_rate_0331_lag, 
                             data = df_ger_death_scaled)

summary(lm_deaths_0331_ctrl)
## 
## Call:
## lm(formula = death_rate_0331 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + death_rate_0331_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8867 -0.4742 -0.2131  0.1644  5.1407 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.003276   0.047481  -0.069 0.945030    
## airport_dist         0.063879   0.060650   1.053 0.292924    
## conservative         0.033146   0.073530   0.451 0.652415    
## male                -0.170458   0.069611  -2.449 0.014805 *  
## age                 -0.282595   0.087552  -3.228 0.001360 ** 
## popdens             -0.056660   0.075956  -0.746 0.456170    
## manufact             0.190338   0.055019   3.460 0.000605 ***
## tourism             -0.071972   0.054777  -1.314 0.189699    
## academics           -0.081969   0.072606  -1.129 0.259654    
## medinc               0.061189   0.065733   0.931 0.352534    
## healthcare          -0.095838   0.062240  -1.540 0.124469    
## temp                -0.120116   0.064268  -1.869 0.062423 .  
## death_rate_0331_lag  0.231548   0.087866   2.635 0.008766 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.924 on 366 degrees of freedom
## Multiple R-squared:  0.1733, Adjusted R-squared:  0.1462 
## F-statistic: 6.392 on 12 and 366 DF,  p-value: 2.503e-10
lm_deaths_0930 <- spec_calculate(df = df_ger_death_scaled, 
               y = "death_rate_0930", 
               model = "lm", 
               controls = covariates %>% append('death_rate_0930_lag'),
               all.comb = T)

calc_all_summaries(lm_deaths_0930)
## $pers_o
##        estimate   std.error statistic  p.value    conf.low  conf.high
## mean 0.09176990 0.055963743 1.6313826 0.145225 -0.01827605 0.20181584
## sd   0.02910357 0.002231346 0.4896016 0.137522  0.02703758 0.03164673
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.12375372         0.1073184 0.94465296      7.626281 0.002057778
## sd      0.03576864         0.0336473 0.01767784      2.103161 0.013217159
##        fit_df  fit_logLik    fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.087716 1042.17543 1077.6133    331.22109      371.000000
## sd   1.732262    7.626127   13.09662   11.0166     13.52055        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2592773        1        0            0.2592773
## sd          0   0.4382916        0        0            0.4382916
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate    std.error  statistic   p.value    conf.low  conf.high
## mean -0.02373502 0.0505987501 -0.4672378 0.6381103 -0.12323134 0.07576129
## sd    0.02148019 0.0009822833  0.4226164 0.2372276  0.02203167 0.02109180
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11782829        0.10128121 0.94783532      7.197791 0.009895088
## sd      0.03662981        0.03456711 0.01808525      2.168291 0.058343172
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.359819 1044.71964 1080.15746    333.46091      371.000000
## sd   1.732262    7.745196   13.38859   11.39777     13.84607        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1352539 0.8647461                    0
## sd          0           0 0.3420363 0.3420363                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05256399 0.051225663 1.0284875 0.3487208 -0.04816508 0.15329305
## sd   0.02270246 0.000942625 0.4474427 0.1962162  0.02322695 0.02232001
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11989229        0.10338102 0.94673451      7.358761 0.003794686
## sd      0.03579107        0.03374005 0.01768661      2.135709 0.021709173
##        fit_df  fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.922157 1043.8443 1079.28214    332.68071      371.000000
## sd   1.732262    7.596441   13.1105   11.22513     13.52902        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379  0.04541016        1        0           0.04541016
## sd          0  0.20822742        0        0           0.20822742
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate  std.error statistic   p.value    conf.low  conf.high
## mean 0.01004584 0.05075540 0.1959914 0.7221615 -0.08975852 0.10985019
## sd   0.02082310 0.00116779 0.4083976 0.1840059  0.02043765 0.02144881
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11736593         0.1008111 0.94808125      7.161593  0.01013347
## sd      0.03683933         0.0347719 0.01819243      2.173253  0.06061994
##        fit_df  fit_logLik    fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.457345 1044.91469 1080.3525    333.63568      371.000000
## sd   1.732262    7.789187   13.46958   11.4449     13.92527        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.6926270 0.3073730                    0
## sd          0           0 0.4614616 0.4614616                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.03733565 0.050718234 0.7335759 0.4914366 -0.06239562 0.13706691
## sd   0.02251156 0.001073727 0.4357044 0.2497344  0.02195929 0.02324317
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11858409        0.10205446 0.94741952      7.243279  0.01153319
## sd      0.03739883        0.03533075 0.01848994      2.202537  0.06886308
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.190383 1044.38077 1079.81859    333.17521      371.000000
## sd   1.732262    7.912561   13.70082   11.59235     14.13676        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379           0 0.9448242 0.05517578                    0
## sd          0           0 0.2283509 0.22835092                    0
##      significant_negative
## mean                    0
## sd                      0
plot_all_curves(lm_deaths_0930, 'lm_deaths_0930')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_deaths_0930_ctrl <- lm(death_rate_0930 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + death_rate_0930_lag, 
                             data = df_ger_death_scaled)

summary(lm_deaths_0930_ctrl)
## 
## Call:
## lm(formula = death_rate_0930 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + death_rate_0930_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9553 -0.5092 -0.2248  0.2556  6.5021 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         -0.004634   0.047233  -0.098  0.92190   
## airport_dist         0.026542   0.060328   0.440  0.66023   
## conservative         0.026796   0.073079   0.367  0.71407   
## male                -0.087034   0.069386  -1.254  0.21052   
## age                 -0.211738   0.086909  -2.436  0.01531 * 
## popdens              0.029944   0.075712   0.396  0.69270   
## manufact             0.135993   0.054799   2.482  0.01353 * 
## tourism             -0.093115   0.054578  -1.706  0.08884 . 
## academics           -0.162248   0.072481  -2.238  0.02579 * 
## medinc               0.093815   0.065408   1.434  0.15234   
## healthcare          -0.026707   0.061800  -0.432  0.66588   
## temp                -0.184670   0.064375  -2.869  0.00436 **
## death_rate_0930_lag  0.284429   0.087564   3.248  0.00127 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9191 on 366 degrees of freedom
## Multiple R-squared:  0.182,  Adjusted R-squared:  0.1552 
## F-statistic: 6.788 on 12 and 366 DF,  p-value: 4.504e-11

Predict Socdist Onset

cox_onset_socdist <- spec_calculate(df = df_ger_socdist_scaled, 
               y = "Surv(cpt_day_socdist_2, event)", 
               model = "cox_model", 
               controls = covariates %>% 
                 append(c('cpt_day_socdist_2_lag', 'onset_prev', 'slope_prev')),
               all.comb = T)

cox_onset_socdist <- cox_onset_socdist %>% map(filter_socdist)

cox_onset_socdist_hr <- cox_onset_socdist %>% map(coef_to_hr)

calc_all_summaries(cox_onset_socdist_hr)
## $pers_o
##        estimate   std.error  statistic   p.value   conf.low  conf.high fit_n
## mean 0.92412115 0.060146732 -1.3205350 0.2201715 0.82135075 1.03977081   379
## sd   0.02156262 0.002248289  0.4019257 0.1387493 0.01899451 0.02532624     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.32341     0.005706523         36.55378
## sd            0          13.75696     0.024395795         12.35949
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.005739844           36.14392      0.005735019                   NA
## sd      0.023985584           12.05654      0.023777053                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09557573         0.9999496      0.61749673
## sd                   NA    0.03296967         0.0000000      0.03336312
##      fit_std.error.concordance   fit_logLik    fit_AIC  fit_BIC fit_nobs
## mean               0.028868295 -1856.052442 3730.10488 3765.543      379
## sd                 0.001210389     6.878482   12.02107   11.250        0
##      significant positive negative significant_positive significant_negative
## mean  0.06811523        0        1                    0           0.06811523
## sd    0.25197430        0        0                    0           0.25197430
## 
## $pers_c
##        estimate   std.error  statistic   p.value  conf.low  conf.high fit_n
## mean 0.92379825 0.053425612 -1.4965299 0.1765796 0.8319259 1.02582137   379
## sd   0.02216087 0.001129165  0.4722299 0.1393372 0.0186606 0.02630592     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.89125     0.002080154         37.09716
## sd            0          12.54155     0.007141225         11.33829
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.002219141           36.66871      0.002244855                   NA
## sd      0.007425369           11.03882      0.007438921                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09703073         0.9999496      0.62825725
## sd                   NA    0.02998817         0.0000000      0.02806221
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean                0.02888723 -1855.768524 3729.53705 3764.97487      379
## sd                  0.00105630     6.270773   10.78844   10.29674        0
##      significant positive negative significant_positive significant_negative
## mean   0.1914062        0        1                    0            0.1914062
## sd     0.3934561        0        0                    0            0.3934561
## 
## $pers_e
##        estimate    std.error  statistic    p.value   conf.low  conf.high fit_n
## mean 0.88750828 0.0562011720 -2.1242714 0.04010842 0.79494811 0.99084875   379
## sd   0.01407696 0.0008904741  0.2702866 0.02493297 0.01325974 0.01506055     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          41.01862     0.002206119         39.25972
## sd            0          13.53017     0.009699107         12.25642
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean     0.00234413           38.82729      0.002322961                   NA
## sd       0.01021985           11.92044      0.010047201                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.1020040         0.9999496      0.62851368
## sd                   NA     0.0321979         0.0000000      0.02755481
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028252338 -1854.704839 3727.40968 3762.84750      379
## sd                 0.001407939     6.765084   11.72497   10.85692        0
##      significant positive negative significant_positive significant_negative
## mean   0.6994629        0        1                    0            0.6994629
## sd     0.4585476        0        0                    0            0.4585476
## 
## $pers_a
##        estimate   std.error  statistic   p.value   conf.low conf.high fit_n
## mean 0.94998162 0.055857806 -0.9252410 0.3798245 0.85145243 1.0599184   379
## sd   0.01686576 0.001206494  0.3290921 0.1672532 0.01424408 0.0201586     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          37.40628      0.00487520         35.73613
## sd            0          13.19216      0.01728311         11.95547
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.005074892           35.31908      0.005112128                   NA
## sd      0.017780079           11.64557      0.017763039                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09343289         0.9999496      0.62513121
## sd                   NA    0.03167007         0.0000000      0.03199302
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028872705 -1856.51101 3731.02201 3766.45984      379
## sd                 0.001126387     6.59608   11.46145   10.84048        0
##      significant positive negative significant_positive significant_negative
## mean           0        0        1                    0                    0
## sd             0        0        0                    0                    0
## 
## $pers_n
##        estimate   std.error statistic    p.value   conf.low  conf.high fit_n
## mean 1.08916744 0.057407948 1.4897329 0.14528486 0.97327553 1.21886580   379
## sd   0.01220407 0.001200459 0.2138959 0.06005665 0.01235991 0.01223976     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.68176      0.00353503         37.02830
## sd            0          13.43207      0.01343931         12.20186
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.003639324           36.63856      0.003605216                   NA
## sd      0.013488312           11.86728      0.013282388                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09645849         0.9999496      0.62898123
## sd                   NA    0.03214506         0.0000000      0.03045481
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028670041 -1855.873268 3729.74654 3765.18436      379
## sd                 0.001247607     6.716033   11.72373   11.09026        0
##      significant positive negative significant_positive significant_negative
## mean           0        1        0                    0                    0
## sd             0        0        0                    0                    0
plot_all_curves(cox_onset_socdist_hr, 'cox_onset_socdist_hr', hr = T)
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

cox_onset_socdist_ctrl <- coxph(Surv(cpt_day_socdist_2, event) ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev + slope_prev +
                                 cpt_day_socdist_2_lag,
                             data = df_ger_socdist_scaled)

summary(cox_onset_socdist_ctrl)
## Call:
## coxph(formula = Surv(cpt_day_socdist_2, event) ~ airport_dist + 
##     conservative + male + age + popdens + manufact + tourism + 
##     academics + medinc + healthcare + temp + onset_prev + slope_prev + 
##     cpt_day_socdist_2_lag, data = df_ger_socdist_scaled)
## 
##   n= 379, number of events= 379 
## 
##                            coef exp(coef)  se(coef)      z Pr(>|z|)    
## airport_dist           0.125877  1.134142  0.067907  1.854   0.0638 .  
## conservative          -0.331611  0.717767  0.082851 -4.002 6.27e-05 ***
## male                  -0.141864  0.867740  0.081693 -1.737   0.0825 .  
## age                   -0.017427  0.982724  0.103089 -0.169   0.8658    
## popdens                0.108749  1.114883  0.091485  1.189   0.2346    
## manufact              -0.072176  0.930367  0.052325 -1.379   0.1678    
## tourism                0.009962  1.010012  0.061504  0.162   0.8713    
## academics             -0.058010  0.943640  0.079147 -0.733   0.4636    
## medinc                -0.104273  0.900979  0.071349 -1.461   0.1439    
## healthcare            -0.173877  0.840400  0.069290 -2.509   0.0121 *  
## temp                   0.077350  1.080420  0.072919  1.061   0.2888    
## onset_prev            -0.013565  0.986527  0.075247 -0.180   0.8569    
## slope_prev             0.146099  1.157311  0.086200  1.695   0.0901 .  
## cpt_day_socdist_2_lag -0.238041  0.788170  0.114252 -2.083   0.0372 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## airport_dist             1.1341     0.8817    0.9928    1.2956
## conservative             0.7178     1.3932    0.6102    0.8443
## male                     0.8677     1.1524    0.7394    1.0184
## age                      0.9827     1.0176    0.8029    1.2028
## popdens                  1.1149     0.8970    0.9319    1.3338
## manufact                 0.9304     1.0748    0.8397    1.0308
## tourism                  1.0100     0.9901    0.8953    1.1394
## academics                0.9436     1.0597    0.8080    1.1020
## medinc                   0.9010     1.1099    0.7834    1.0362
## healthcare               0.8404     1.1899    0.7337    0.9626
## temp                     1.0804     0.9256    0.9365    1.2464
## onset_prev               0.9865     1.0137    0.8513    1.1433
## slope_prev               1.1573     0.8641    0.9774    1.3703
## cpt_day_socdist_2_lag    0.7882     1.2688    0.6300    0.9860
## 
## Concordance= 0.678  (se = 0.029 )
## Likelihood ratio test= 61.37  on 14 df,   p=7e-08
## Wald test            = 57.4  on 14 df,   p=3e-07
## Score (logrank) test = 58.42  on 14 df,   p=2e-07

Predict Socdist Mean

lm_mean_socdist <- spec_calculate(df = df_ger_socdist_scaled, 
               y = "mean_diff_socdist_2", 
               model = "lm", 
               controls = covariates %>% 
                 append(c('mean_diff_socdist_2_lag', 'onset_prev', 'slope_prev')), 
               all.comb = T)

lm_mean_socdist <- lm_mean_socdist %>% map(filter_socdist)

calc_all_summaries(lm_mean_socdist)
## $pers_o
##        estimate   std.error statistic    p.value   conf.low  conf.high
## mean 0.15330885 0.044355910  3.425476 0.02297009 0.06608679 0.24053090
## sd   0.07755911 0.001414087  1.656772 0.03632131 0.07585609 0.07932302
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean    0.45206387        0.43888864 0.7476768     35.209876 7.055131e-20
## sd      0.07030338        0.07031662 0.0457296      7.949236 2.044565e-18
##        fit_df fit_logLik   fit_AIC  fit_BIC fit_deviance fit_df.residual
## mean 9.000000  -421.8148 865.62962 908.9425    207.11986      369.000000
## sd   1.732262    23.1611  44.31669  40.9420     26.57468        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.8122559        1        0            0.8122559
## sd          0   0.3905554        0        0            0.3905554
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate  std.error  statistic   p.value    conf.low  conf.high
## mean -0.03851258 0.04069143 -0.9421542 0.4218812 -0.11852873 0.04150356
## sd    0.02790481 0.00308532  0.6560585 0.2885431  0.02975273 0.02730813
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.43015918        0.41651193 0.76145621     32.660794 8.760185e-07
## sd      0.09616314        0.09669053 0.06060868      9.707367 2.157079e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.28432 878.56865 921.88154    215.39983      369.000000
## sd   1.732262   29.50691  57.06175  53.66655     36.34967        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.0769043 0.0546875 0.9453125                    0
## sd          0   0.2664721 0.2273970 0.2273970                    0
##      significant_negative
## mean            0.0769043
## sd              0.2664721
## 
## $pers_e
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.01503002 0.041886341 -0.3896659 0.5811354 -0.09739585 0.06733580
## sd    0.02278655 0.003250974  0.5237195 0.1998854  0.01851793 0.02787933
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.42867767        0.41499876 0.76240500      32.47312 3.773528e-06
## sd      0.09720623        0.09774027 0.06116183       9.72353 8.896818e-05
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.73869 879.4774 922.79028    215.95984      369.000000
## sd   1.732262   29.73164  57.4969  54.06717     36.74395        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1518555 0.8481445                    0
## sd          0           0 0.3589246 0.3589246                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.02332656 0.040841977 -0.5897552 0.4122219 -0.10363874 0.05698563
## sd    0.04379788 0.002890095  1.0885441 0.2944594  0.04279858 0.04549033
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.43026173        0.41662073 0.76130558     32.712822 6.230705e-06
## sd      0.09788511        0.09844834 0.06159559      9.858796 1.500141e-04
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.17439 878.3488 921.66168    215.36107      369.000000
## sd   1.732262   29.92968  57.9119  54.51971     37.00057        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1582031 0.3818359 0.6181641                    0
## sd          0   0.3649759 0.4858960 0.4858960                    0
##      significant_negative
## mean            0.1582031
## sd              0.3649759
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.00604797 0.041254131 -0.1365582 0.7023933 -0.08717062 0.07507468
## sd    0.01968351 0.003240101  0.4720721 0.2018550  0.02250371 0.01869854
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.42841045        0.41472473 0.76258917     32.435635 6.023785e-06
## sd      0.09718258        0.09771424 0.06110589      9.710442 1.448236e-04
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.83354 879.66708 922.97998    216.06085      369.000000
## sd   1.732262   29.68545  57.40582  53.97958     36.73502        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.4108887 0.5891113                    0
## sd          0           0 0.4920552 0.4920552                    0
##      significant_negative
## mean                    0
## sd                      0
plot_all_curves(lm_mean_socdist, 'lm_mean_socdist')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_mean_socdist_ctrl <- lm(mean_diff_socdist_2 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev + slope_prev +
                                 mean_diff_socdist_2_lag,
                             data = df_ger_socdist_scaled)

summary(lm_mean_socdist_ctrl)
## 
## Call:
## lm(formula = mean_diff_socdist_2 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + onset_prev + slope_prev + mean_diff_socdist_2_lag, 
##     data = df_ger_socdist_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.85241 -0.44592  0.01798  0.37505  2.81498 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -0.001481   0.035208  -0.042   0.9665    
## airport_dist            -0.010038   0.045187  -0.222   0.8243    
## conservative            -0.273373   0.055239  -4.949 1.14e-06 ***
## male                    -0.117690   0.051800  -2.272   0.0237 *  
## age                      0.003341   0.067934   0.049   0.9608    
## popdens                  0.328194   0.056804   5.778 1.63e-08 ***
## manufact                 0.085452   0.041025   2.083   0.0380 *  
## tourism                 -0.013999   0.041163  -0.340   0.7340    
## academics                0.364853   0.054490   6.696 8.15e-11 ***
## medinc                  -0.073709   0.050032  -1.473   0.1416    
## healthcare              -0.008339   0.046183  -0.181   0.8568    
## temp                    -0.081000   0.048476  -1.671   0.0956 .  
## onset_prev               0.008089   0.051930   0.156   0.8763    
## slope_prev               0.031435   0.055945   0.562   0.5745    
## mean_diff_socdist_2_lag  0.135695   0.064960   2.089   0.0374 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6853 on 364 degrees of freedom
## Multiple R-squared:  0.5478, Adjusted R-squared:  0.5304 
## F-statistic: 31.49 on 14 and 364 DF,  p-value: < 2.2e-16

Save all results

results_ger <- list(cox_onset_prev_hr, lm_slope_prev, lm_slope_prev_var, 
                   lm_cases_0331, lm_cases_0930, lm_deaths_0331, lm_deaths_0930, 
                   cox_onset_socdist_hr, lm_mean_socdist)


names(results_ger) <- list('cox_onset_prev_hr', 'lm_slope_prev', 'lm_slope_prev_var', 
                   'lm_cases_0331', 'lm_cases_0930', 'lm_deaths_0331', 'lm_deaths_0930', 
                   'cox_onset_socdist_hr', 'lm_mean_socdist')

save(results_ger, file="GER_spec_results.RData")

Descriptives

Distributions

df_ger_desc <- df_ger_socdist %>% 
  dplyr::select(pers_o, pers_c, pers_e, pers_a, pers_n, 
                age, male, conservative, 
                academics, medinc, manufact, 
                airport_dist, tourism, healthcare, popdens) %>% 
  distinct() %>% drop_na()

ger_means <- df_ger_desc %>% summarise_all(mean)
ger_sd <- df_ger_desc %>% summarise_all(sd)

desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame()

desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger  

Explore correlations

a <- df_ger_socdist_unscaled %>% select(kreis, pers_o:pers_n, cpt_day_socdist, mean_diff_socdist)
b <- df_ger_prev_unscaled %>% select(kreis, onset_prev, slope_prev)
c <- df_ger_death_unscaled %>% select(kreis, case_rate_0331:death_rate_0930)

ger_joined <- plyr::join_all(list(a, b, c), by='kreis')

ger_joined %>% select(-kreis) %>% cor(use = 'pairwise.complete')
##                        pers_o      pers_c      pers_e      pers_a        pers_n
## pers_o             1.00000000 -0.10201270  0.27983871  0.15759347 -0.0640347187
## pers_c            -0.10201270  1.00000000  0.21472242  0.32690021 -0.3551086811
## pers_e             0.27983871  0.21472242  1.00000000  0.18288787 -0.3649780097
## pers_a             0.15759347  0.32690021  0.18288787  1.00000000 -0.4280315714
## pers_n            -0.06403472 -0.35510868 -0.36497801 -0.42803157  1.0000000000
## cpt_day_socdist   -0.10416238  0.10437674 -0.04952867 -0.01672102  0.0138442293
## mean_diff_socdist  0.31429946 -0.12682312  0.09840937 -0.08341683  0.0006497013
## onset_prev        -0.16574985  0.05702967 -0.27615283 -0.09470167  0.2352513229
## slope_prev         0.15470785 -0.04204164  0.24306903  0.13157797 -0.1983834025
## case_rate_0331     0.10695243 -0.01752091  0.20572870  0.13949525 -0.1722746724
## death_rate_0331    0.06566400 -0.03556727  0.10805800  0.09283184 -0.0999658180
## case_rate_0930     0.15634007 -0.05283919  0.25572247  0.04675307 -0.1116940920
## death_rate_0930    0.07587651 -0.03740828  0.08958204  0.02599047 -0.0095583322
##                   cpt_day_socdist mean_diff_socdist  onset_prev  slope_prev
## pers_o              -0.1041623824      0.3142994608 -0.16574985  0.15470785
## pers_c               0.1043767434     -0.1268231184  0.05702967 -0.04204164
## pers_e              -0.0495286726      0.0984093720 -0.27615283  0.24306903
## pers_a              -0.0167210202     -0.0834168254 -0.09470167  0.13157797
## pers_n               0.0138442293      0.0006497013  0.23525132 -0.19838340
## cpt_day_socdist      1.0000000000      0.0358262655  0.14392128 -0.11996515
## mean_diff_socdist    0.0358262655      1.0000000000 -0.20519894  0.22948985
## onset_prev           0.1439212762     -0.2051989385  1.00000000 -0.67617865
## slope_prev          -0.1199651501      0.2294898547 -0.67617865  1.00000000
## case_rate_0331      -0.0689721829      0.1733238658 -0.66024501  0.92767374
## death_rate_0331      0.0003684063      0.1435286041 -0.40795429  0.69587384
## case_rate_0930      -0.1260015001      0.2185131719 -0.58477976  0.82060844
## death_rate_0930     -0.0078242968      0.1456354879 -0.32697005  0.64249969
##                   case_rate_0331 death_rate_0331 case_rate_0930 death_rate_0930
## pers_o                0.10695243    0.0656640020     0.15634007     0.075876508
## pers_c               -0.01752091   -0.0355672682    -0.05283919    -0.037408282
## pers_e                0.20572870    0.1080579960     0.25572247     0.089582044
## pers_a                0.13949525    0.0928318364     0.04675307     0.025990472
## pers_n               -0.17227467   -0.0999658180    -0.11169409    -0.009558332
## cpt_day_socdist      -0.06897218    0.0003684063    -0.12600150    -0.007824297
## mean_diff_socdist     0.17332387    0.1435286041     0.21851317     0.145635488
## onset_prev           -0.66024501   -0.4079542869    -0.58477976    -0.326970047
## slope_prev            0.92767374    0.6958738421     0.82060844     0.642499692
## case_rate_0331        1.00000000    0.7967609328     0.82553862     0.687015730
## death_rate_0331       0.79676093    1.0000000000     0.68600498     0.850162063
## case_rate_0930        0.82553862    0.6860049833     1.00000000     0.745023167
## death_rate_0930       0.68701573    0.8501620632     0.74502317     1.000000000
ger_means <- ger_joined %>% select(-kreis) %>% summarise_all(mean)
ger_sd <- ger_joined %>% select(-kreis) %>% summarise_all(sd)

desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame() %>%
  rename(mean=V1, sd=V2)

desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger  

Visualizations

df_ger_prev <- read_csv('GER_prevalence.csv')

df_ger_prev <- df_ger_prev %>%
  mutate(date = as.Date(date, "%d%b%Y"),
         kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day) %>%
  filter(date < '2020-04-01')

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')

tto <- df_ger_prev_unscaled %>%
  filter(event==1) %>%
  merge(df_dates, by.x = 'onset_prev', by.y = 'time') %>%
  ggplot(aes(x=date+1)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - COVID-19 Onset') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Date')
  
soa <- df_ger_slope_prev %>% 
  ggplot(aes(x=slope_prev)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - COVID-19 Growth Rates') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Standardized growth rates')

figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)

figure

tto <- df_ger_death %>% 
  ggplot(aes(x=case_rate_0930)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Cumulative Case Rates') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Cases Per 1000 Inhabitants')
  
soa <- df_ger_death %>%
  ggplot(aes(x=death_rate_0930)) +
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Cumulative Death Rates') +
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Deaths Per 1000 Inhabitants')


figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)
figure

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')


tto <- df_ger_cpt_socdist %>% 
  merge(df_dates, by.x = 'cpt_day_socdist_2', by.y = 'time') %>%
  #filter(cpt_day_socdist_2 < 30) %>% 
  ggplot(aes(x=date+1)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Time to Adoption') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Date') +
  xlim(c(as.Date('2020-03-01'), as.Date('2020-03-31')))
  
soa <- df_ger_cpt_socdist %>% 
  ggplot(aes(x=mean_diff_socdist_2)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Strength of Adjustment') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Standardized mean difference')


figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)
figure

======= COVID19 GER

Prepare data

Read and format prevalence data

df_ger_prev <- read_csv('GER_prevalence.csv')

df_ger_prev <- df_ger_prev %>%
  mutate(date = as.Date(date, "%d%b%Y"),
         kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day) %>%
  filter(date < '2020-04-01')

df_ger_prev %>% write_csv('ger_prev_kreis.csv')
df_ger_death <- read_csv('GER_prevalence.csv')

df_ger_death <-  df_ger_death %>%
  mutate(date = as.Date(date, "%d%b%Y"), 
         kreis = as.character(kreis))

df_ger_death_0930 <- df_ger_death %>% 
  filter(date == '2020-09-30') %>% 
  rename(case_rate_0930 = rate_day,
         death_rate_0930 = death_day) %>% 
  select(kreis, case_rate_0930, death_rate_0930)

df_ger_death_0331 <- df_ger_death %>% 
  filter(date == '2020-03-31') %>% 
  rename(case_rate_0331 = rate_day,
         death_rate_0331 = death_day) %>% 
  select(kreis, case_rate_0331, death_rate_0331)


df_ger_death <- merge(df_ger_death_0331, df_ger_death_0930)

# save df
df_ger_death %>% write_csv('ger_death_kreis_maps.csv')

Read and format scoial distancing data

fb_files <- list.files('../FB Data/GER individual files',
                       '*.csv', full.names = T)

df_ger_socdist <- fb_files %>% 
  map(read.csv) %>% bind_rows()

kreis_match <- read_csv('kreismatch.csv')

df_ger_socdist <- df_ger_socdist %>%
  select(-polygon_name) %>%
  plyr::join(kreis_match, by = 'polygon_id') %>%
  select(kreis, ds, all_day_bing_tiles_visited_relative_change,
         all_day_ratio_single_tile_users) %>%
  rename(date = ds,
         socdist_tiles = all_day_bing_tiles_visited_relative_change,
         socdist_single_tile = all_day_ratio_single_tile_users) %>%
  mutate(kreis = as.character(kreis),
         date = as.Date(date)) %>%
  arrange(kreis, date)

Read and format personality data

df_ger_pers <- read_csv('GER_personality_2.csv')

df_ger_pers %>% filter(frequ >= 50) %>% .$frequ %>% mean()
## [1] 284.0416
df_ger_pers <- df_ger_pers %>% 
  filter(frequ >= 50) %>%
  select(kreis, open, sci, extra, agree, neuro) %>%
  dplyr::rename(pers_o = open,
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct() %>%
  mutate(kreis = as.character(kreis))

Read and format county level controls

df_ger_ctrl <- read_csv('GER_controls_new.csv')

df_ger_ctrl <- df_ger_ctrl %>%
  select(-gdp, -cdu, -kreis_nme) %>%
  rename(tourism = tourism_beds,
         healthcare = hospital_beds,
         airport_dist = airport,
         conservative = afd, 
         popdens = popdens_new) %>%
  mutate(kreis = as.character(kreis))


df_ger_temp <- read_csv("GER_controls_weather_2.csv") %>% 
  select(kreis, kr_tamm_202003) %>% 
  rename(temp = kr_tamm_202003) %>%
  mutate(kreis = as.character(as.numeric(kreis)))

df_ger_ctrl <- df_ger_ctrl %>% left_join(df_ger_temp, 'kreis')

df_ger_pers %>% merge(df_ger_ctrl, by='kreis') %>% 
  write_csv('df_ger_pers_kreis.csv')

Merge prevalence data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge prevalence data
df_ger_prev <- df_ger_prev %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>%
  merge(df_dates, by='date') %>%
  arrange(kreis, date) %>%
  drop_na()

df_ger_prev %>% select(kreis) %>% distinct() %>% nrow()
## [1] 383
# join data frames
df_ger_death <- df_ger_death %>%
  plyr::join(df_ger_ctrl, by='kreis') %>%
  plyr::join(df_ger_pers, by='kreis') %>%
  drop_na()

# save df
df_ger_death %>%
  #select(kreis, case_rate, death_rate) %>%
  write_csv('ger_death_kreis_controls.csv')

Merge social distancing data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')

# merge socdist data
df_ger_socdist <- df_ger_socdist %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>% 
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_socdist %>% select(kreis) %>% distinct() %>% nrow()
## [1] 379

Control for weekend effect

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)

df_ger_loess <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_loess_2 <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_tiles ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  rename(loess_2 = loess) %>%
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_socdist <- df_ger_socdist %>% 
  merge(df_ger_loess, by=c('kreis', 'time')) %>% 
  merge(df_ger_loess_2, by=c('kreis', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile),
         socdist_tiles_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess_2, socdist_tiles)) %>%
  arrange(kreis, time) %>% 
  select(-weekday)

df_ger_socdist <- df_ger_socdist %>% drop_na() %>% mutate(time = time-1)

Include cleaned social distancing data

df_ger_socdist <- df_ger_socdist %>% 
  mutate(socdist_single_tile = socdist_single_tile_clean,
         socdist_tiles = socdist_tiles_clean) %>% 
  select(-loess, -loess_2, -socdist_single_tile_clean, -socdist_tiles_clean)

Define outcomes of interest

COVID - Extract onset

# get onset day
df_ger_onset_prev <- df_ger_prev %>% 
  group_by(kreis) %>% 
  filter(rate_day > 0.1) %>%
  summarise(onset_prev = min(time))
  
# merge with county data
df_ger_onset_prev <- df_ger_prev %>% 
  select(-date, -time, -rate_day) %>%
  distinct() %>% 
  left_join(df_ger_onset_prev, by = 'kreis')

# handle censored data
df_ger_onset_prev <- df_ger_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_ger_prev$date)))+1))

COVID - Extract slopes

# cut time series
df_ger_prev <- df_ger_prev %>% 
  filter(date > '2020-03-01' & date < '2020-04-01')

# extract slope prevalence
df_ger_slope_prev <- df_ger_prev %>% 
  split(.$kreis) %>% 
  map(~ lm(log(rate_day+1) ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>%
  rename(slope_prev = '.')

# merge with county data
df_ger_slope_prev <- df_ger_onset_prev %>% 
  left_join(df_ger_slope_prev, by = 'kreis')

# get unscaled object for descriptives
df_ger_prev_desc <- df_ger_slope_prev
# cut time series before onset
df_ger_slope_var <- df_ger_prev %>% 
  filter(date <= '2020-04-30') %>%
  group_by(kreis) %>% 
  filter(rate_day > 0.1) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30)

# extract slope prevalence
df_ger_slope_var <- df_ger_slope_var %>% 
  split(.$kreis) %>% 
  map(~ lm(log(rate_day+1) ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>% 
  rename(slope_prev_var = '.')

# merge with county data
df_ger_slope_prev <- df_ger_slope_prev %>% 
  left_join(df_ger_slope_var, by = 'kreis') %>%
  mutate(slope_prev_var = replace_na(slope_prev_var, 0))

Social Distancing - Change point analysis

# keep only counties with full data
kreis_complete <- df_ger_socdist %>% 
  group_by(kreis) %>% 
  summarise(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$kreis

# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist %>% 
  select(kreis, socdist_single_tile) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

df_ger_socdist_cpt_results_2 <- df_ger_socdist %>% 
  select(kreis, socdist_tiles) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_tiles),
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('kreis')

df_ger_socdist_cpt_day_2 <- df_ger_socdist_cpt_results_2 %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist_2 = '.') %>%
  rownames_to_column('kreis')

# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('kreis')

df_ger_socdist_cpt_mean_diff_2 <- df_ger_socdist_cpt_results_2 %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ -.[2]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist_2 = '.') %>%
  rownames_to_column('kreis')

# calculate means 
df_ger_socdist_mean <- df_ger_socdist %>%
  group_by(kreis) %>%
  summarise(mean_socdist = mean(socdist_single_tile))

df_ger_socdist_mean_2 <- df_ger_socdist %>%
  group_by(kreis) %>%
  summarise(mean_socdist_2 = -mean(socdist_tiles))

# merge with county data
df_ger_cpt_socdist <- df_ger_socdist %>%
  select(-date, -time, -socdist_single_tile, -socdist_tiles) %>%
  distinct() %>%
  left_join(df_ger_socdist_cpt_day, by='kreis') %>%
  left_join(df_ger_socdist_cpt_day_2, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff_2, by='kreis') %>%
  left_join(df_ger_socdist_mean, by='kreis') %>%
  left_join(df_ger_socdist_mean_2, by='kreis') %>%
  left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
  left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis')

# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), 
                                  as.numeric(diff(range(df_ger_socdist$date))), 
                                  cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 
                          as.numeric(diff(range(df_ger_socdist$date))), 0, 1))

Remove incomplete cases

df_ger_prev_unscaled <- inner_join(df_ger_slope_prev, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_death_unscaled <- inner_join(df_ger_death, df_ger_cpt_socdist %>% select(kreis), by='kreis')
df_ger_socdist_unscaled <- inner_join(df_ger_cpt_socdist, df_ger_slope_prev %>% select(kreis), by='kreis')

kreis_list <- df_ger_prev_unscaled %>% select(kreis)

kreis_list %>% write_csv('kreis_list.csv')

df_ger_prev_unscaled %>% merge(df_ger_death_unscaled %>% select(kreis:death_rate_0930))  %>% 
  merge(df_ger_socdist_unscaled %>% select(kreis, cpt_day_socdist_2, mean_diff_socdist_2)) %>% 
  select(-event) %>% write_csv('ger_all_variables.csv')

Rescale Data

df_ger_prev_scaled <- df_ger_prev_unscaled %>% 
  mutate_at(vars(-kreis, -event), scale)
df_ger_death_scaled <- df_ger_death_unscaled %>% 
  mutate_at(vars(-kreis), scale)
df_ger_socdist_scaled <- df_ger_socdist_unscaled %>%
  mutate_at(vars(-kreis, -event), scale)

Calculate spatial lags

# create UDF to calculate lagged variables
calc_lags <- function(df, weights, cols){
  
  cols_only <- df %>% select(cols)
  cols_lag <- weights %*% as.matrix(cols_only) %>% 
  as.matrix() %>% as.data.frame()
names(cols_lag) <- paste0(names(cols_lag), '_lag')

return(cols_lag)
}
# read_weight matrix 
weight_mat_norm <- read_csv('GER_spatial_weights.csv')

weight_mat_norm <- weight_mat_norm %>% select(-kreis) %>% as.matrix()

dim(weight_mat_norm)
## [1] 379 379
# generate spatially lagged y 
y_only <- df_ger_prev_scaled %>% select(onset_prev,slope_prev,slope_prev_var) %>% names()
y_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_prev_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_prev_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_prev_scaled <- cbind(df_ger_prev_scaled, y_lag, X_lag)
# generate spatially lagged y 
y_only <- df_ger_death_scaled %>% select(case_rate_0331:death_rate_0930) %>% names()
y_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_death_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_death_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_death_scaled <- cbind(df_ger_death_scaled, y_lag, X_lag)
# generate spatially lagged y 
y_only <- df_ger_socdist_scaled %>% select(cpt_day_socdist_2,mean_diff_socdist_2) %>% names()
y_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, y_only)

# generate spatially lagged X
X_only <- df_ger_socdist_scaled %>% select(academics:pers_n) %>% names()
X_lag <- calc_lags(df_ger_socdist_scaled, weight_mat_norm, X_only)

# bind new variables to df
df_ger_socdist_scaled <- cbind(df_ger_socdist_scaled, y_lag, X_lag)
write_csv(df_ger_prev_scaled, 'df_ger_slope_prev.csv')
write_csv(df_ger_death_scaled, 'df_ger_death_scaled.csv')
write_csv(df_ger_socdist_scaled, 'df_ger_cpt_socdist.csv')

Run Specification Curve Analysis

# define function to calculate specification curve analyses for all traits
spec_calculate <- function(df, y, model, controls, all.comb = T){
  
  spec_results_o <- run_specs(df = df,
                       y = c(y), x = c("pers_o"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_c <- run_specs(df = df, 
                       y = c(y), x = c("pers_c"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_e <- run_specs(df = df, 
                       y = c(y), x = c("pers_e"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_a <- run_specs(df = df, 
                       y = c(y), x = c("pers_a"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results_n <- run_specs(df = df, 
                       y = c(y), x = c("pers_n"), 
                       model = c(model), controls = controls, 
                       all.comb = all.comb)
  
  spec_results <- list(spec_results_o, spec_results_c, spec_results_e, 
                      spec_results_a, spec_results_n)
  
  names(spec_results) <- list('spec_results_o', 'spec_results_c', 'spec_results_e', 
                      'spec_results_a', 'spec_results_n')

  return(spec_results)
  }
# adapted from specr package code
format_results <- function(df, var, null = 0, desc = FALSE) {

  # rank specs
  if (isFALSE(desc)) {
    df <- df %>%
      dplyr::arrange(!! var)
  } else {
    df <- df %>%
      dplyr::arrange(desc(!! var))
  }

  # create rank variable and color significance
  df <- df %>%
    dplyr::mutate(specifications = 1:nrow(df),
                  color = case_when(conf.low > null ~ "#377eb8",
                                    conf.high < null ~ "#e41a1c",
                                    TRUE ~ "darkgrey"))
  return(df)
}


# define function to plot single specification curve
plot_curves <- function(results, filename, hr=F){
  
  file_path_eps <- paste0("../../Plots/GER/", filename,".eps")
  file_path_pdf <- paste0("../../Plots/GER/", filename,".pdf")
  file_path_png <- paste0("../../Plots/GER/", filename,".png")

  if(hr==F){
    
    results %>%
    format_results(var = .$estimate, null = 0, desc = F) %>%
    ggplot(aes(x = specifications,
               y = estimate,
               ymin = conf.low,
               ymax = conf.high,
               color = color)) +
    geom_point(aes(color = color),
               size = 1) +
    theme_minimal() +
    scale_color_identity() +
    theme(strip.text = element_blank(),
          axis.line = element_line("black", size = .5),
          legend.position = "none",
          panel.spacing = unit(.75, "lines"),
          axis.text = element_text(colour = "black")) +
    labs(x = "") +
      geom_pointrange(alpha = 0.05,
                      size = .6,
                      fatten = 1) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    labs(x = "", y = "standarized coefficient") + 
      coord_fixed(ratio = 2000, ylim = c(-0.4, 0.4)) +
      theme(panel.border = element_blank(), panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            axis.line = element_line(colour = "black")) +
      ggsave(file=file_path_eps, device = 'eps')+
      ggsave(file=file_path_pdf, device = 'pdf')+
      ggsave(file=file_path_png, device = 'png')
    
  }else{
    
    results %>%
    format_results(var = .$estimate, null = 1, desc = F) %>%
    ggplot(aes(x = specifications,
               y = estimate,
               ymin = conf.low,
               ymax = conf.high,
               color = color)) +
    geom_point(aes(color = color),
               size = 1) +
    theme_minimal() +
    scale_color_identity() +
    theme(strip.text = element_blank(),
          axis.line = element_line("black", size = .5),
          legend.position = "none",
          panel.spacing = unit(.75, "lines"),
          axis.text = element_text(colour = "black")) +
    labs(x = "") +
      geom_pointrange(alpha = 0.05,
                      size = .6,
                      fatten = 1) +
    geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
    labs(x = "", y = "standarized coefficient") + 
      coord_fixed(ratio = 2000, ylim = c(0.6, 1.4)) +
      theme(panel.border = element_blank(), panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(), 
      axis.line = element_line(colour = "black")) +
      ggsave(file=file_path_eps, device = 'eps')+
      ggsave(file=file_path_pdf, device = 'pdf')+
      ggsave(file=file_path_png, device = 'png')  
    }
  
}

# define function to plot all curves in list
plot_all_curves <- function(ls, analysis, hr=F){
  
  pers <- c('o', 'c', 'e', 'a', 'n')
  filenames <- as.list(paste0(analysis, '_', pers))
  hr <- as.list(rep(hr, 5))
  pmap(list(ls, filenames, hr), plot_curves)
}
# define function to calculate summary stats for single specification curve
calc_summary <- function(df){
  
  dft <- df %>% select(estimate:fit_nobs)
  dft <- dft %>% mutate(significant = as.numeric(p.value<0.05),
                        positive = as.numeric(statistic>0), 
                        negative = as.numeric(statistic<0), 
                        significant_positive = as.numeric(p.value<0.05 & statistic>0),
                        significant_negative = as.numeric(p.value<0.05 & statistic<0))
  
  mean_temp <- dft %>% map_if(is.numeric, mean, .else=NULL) %>% as.data.frame()
  sd_temp <- dft %>% map_if(is.numeric, sd, .else=NULL) %>% as.data.frame()
  
  df_temp <- rbind(mean_temp, sd_temp)
  row.names(df_temp) <- c('mean', 'sd')
  
  return(df_temp)
}

# define function to calculate summary stats for all curves in list
calc_all_summaries <- function(ls){
  
  ls <- ls %>% map(calc_summary)
  names(ls) <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
  return(ls)
}
coef_to_hr <- function(df){
  
  df %>% mutate(conf.low = exp(estimate-1.96*std.error),
           conf.high = exp(estimate+1.96*std.error),
           estimate = exp(estimate)) 
}

filter_socdist <- function(df){
  
  df %>% filter(str_detect(controls, 'onset_prev') & 
                  str_detect(controls, 'slope_prev'))
}
covariates <- c("airport_dist", "conservative", 'male', 'age', 'popdens',
                'manufact', 'tourism', 'academics', 'medinc', 'healthcare',
                'temp')

Predict Onset

cox_model <- function(formula, data){
  formula <- as.formula(formula)
  coxph(formula = formula, data = data)}
cox_onset_prev <- spec_calculate(df = df_ger_prev_scaled, 
               y = "Surv(onset_prev, event)", 
               model = "cox_model", 
               controls = covariates %>% append('onset_prev_lag'),
               all.comb = T)

cox_onset_prev_hr <- cox_onset_prev %>% map(coef_to_hr)

calc_all_summaries(cox_onset_prev_hr)
## $pers_o
##        estimate   std.error statistic   p.value   conf.low  conf.high fit_n
## mean 1.02726691 0.055320547 0.4779412 0.5231112 0.92180104 1.14483537   379
## sd   0.05390301 0.002840416 0.9535330 0.2800335 0.05018495 0.05845821     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         131.39127    3.811112e-06         140.8145
## sd            0          36.26568    5.698238e-05          41.9226
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   4.156490e-06          137.52997     4.234893e-06                   NA
## sd     5.770949e-05           40.56007     5.819916e-05                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.2895930         0.9999496       0.7195157
## sd                   NA     0.0722395         0.0000000       0.0381976
##      fit_std.error.concordance  fit_logLik   fit_AIC    fit_BIC fit_nobs
## mean               0.015635755 -1809.51851 3633.0370 3660.59978      379
## sd                 0.001317056    18.13284   34.3554   31.38233        0
##      significant  positive  negative significant_positive significant_negative
## mean  0.08911133 0.6496582 0.3503418           0.08911133                    0
## sd    0.28493915 0.4771352 0.4771352           0.28493915                    0
## 
## $pers_c
##        estimate   std.error  statistic   p.value   conf.low  conf.high fit_n
## mean 0.94922252 0.053595445 -0.9954995 0.3479716 0.85456571 1.05437271   379
## sd   0.04503002 0.001467328  0.8933022 0.2887067 0.04054931 0.05018395     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         132.01667    0.0001954177         140.8434
## sd            0          37.85635    0.0029957287          43.1195
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   0.0001852582          137.59893     0.0001881611                   NA
## sd     0.0028016025           41.78484     0.0028121428                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.29044877         0.9999496      0.71826746
## sd                   NA    0.07571748         0.0000000      0.04121107
##      fit_std.error.concordance  fit_logLik    fit_AIC   fit_BIC fit_nobs
## mean               0.015744358 -1809.20581 3632.41162 3659.9744      379
## sd                 0.001364156    18.92817   35.97484   33.0281        0
##      significant  positive  negative significant_positive significant_negative
## mean   0.1445312 0.1584473 0.8415527                    0            0.1445312
## sd     0.3516705 0.3652045 0.3652045                    0            0.3516705
## 
## $pers_e
##        estimate   std.error statistic   p.value   conf.low  conf.high fit_n
## mean 1.10484855 0.052014824 1.9131109 0.1418816 0.99784541 1.22333370   379
## sd   0.05117709 0.001260917 0.9230424 0.1475011 0.04805772 0.05449954     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         134.80506    6.667657e-08        144.44569
## sd            0          33.69771    7.667802e-07         39.58589
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   1.086854e-07          141.00631     1.184551e-07                   NA
## sd     1.098564e-06           38.16651     1.180765e-06                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.29642439         0.9999496      0.72680567
## sd                   NA    0.06635634         0.0000000      0.03214107
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.015317233 -1807.81162 3629.62324 3657.18599      379
## sd                 0.001142571    16.84885   31.78596   28.87968        0
##      significant positive negative significant_positive significant_negative
## mean   0.3879395        1        0            0.3879395                    0
## sd     0.4873401        0        0            0.4873401                    0
## 
## $pers_a
##      estimate   std.error statistic  p.value   conf.low  conf.high fit_n
## mean 1.006651 0.053837899 0.1292402 0.645567 0.90592390 1.11859211   379
## sd   0.032338 0.001749636 0.6071196 0.238925 0.03165915 0.03294262     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          130.6175    0.0006586847        140.37322
## sd            0           37.6386    0.0123791555         43.20648
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   0.0006630851          137.15705      0.000668293                   NA
## sd     0.0122500086           41.79493      0.012237339                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.2878630         0.9999496      0.71879483
## sd                   NA     0.0756924         0.0000000      0.04050933
##      fit_std.error.concordance fit_logLik   fit_AIC    fit_BIC fit_nobs
## mean               0.015616239 -1809.9054 3633.8108 3661.37351      379
## sd                 0.001276486    18.8193   35.7401   32.76038        0
##      significant  positive  negative significant_positive significant_negative
## mean           0 0.4870605 0.5129395                    0                    0
## sd             0 0.4998936 0.4998936                    0                    0
## 
## $pers_n
##       estimate    std.error  statistic    p.value   conf.low  conf.high fit_n
## mean 0.8618104 0.0536108475 -2.7813176 0.01334193 0.77585634 0.95728861   379
## sd   0.0256510 0.0006767025  0.5515346 0.01683804 0.02333785 0.02824371     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379         138.42138    8.297546e-07        148.34185
## sd            0          35.84265    1.236878e-05         41.64845
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean   9.623814e-07          144.57863     1.004892e-06                   NA
## sd     1.340854e-05           40.36034     1.343165e-05                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.30271384         0.9999496      0.72815172
## sd                   NA    0.07048457         0.0000000      0.03346175
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.015109031 -1806.00346 3626.00692 3653.56967      379
## sd                 0.001124059    17.92132   33.95795   31.05319        0
##      significant positive negative significant_positive significant_negative
## mean   0.9492188        0        1                    0            0.9492188
## sd     0.2195775        0        0                    0            0.2195775
plot_all_curves(cox_onset_prev_hr, 'cox_onset_prev_hr', hr = T)
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev_lag, 
                             data = df_ger_prev_scaled)

summary(cox_onset_prev_ctrl)
## Call:
## coxph(formula = Surv(onset_prev, event) ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + onset_prev_lag, data = df_ger_prev_scaled)
## 
##   n= 379, number of events= 379 
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)    
## airport_dist   -0.09147   0.91259  0.07045 -1.298 0.194133    
## conservative   -0.08051   0.92264  0.08006 -1.006 0.314608    
## male           -0.13923   0.87003  0.07451 -1.869 0.061664 .  
## age            -0.25167   0.77750  0.09118 -2.760 0.005779 ** 
## popdens        -0.09156   0.91251  0.07416 -1.235 0.216965    
## manufact        0.00861   1.00865  0.05321  0.162 0.871458    
## tourism        -0.11689   0.88968  0.05988 -1.952 0.050933 .  
## academics       0.07279   1.07551  0.08113  0.897 0.369583    
## medinc          0.29682   1.34558  0.07898  3.758 0.000171 ***
## healthcare     -0.06525   0.93683  0.06938 -0.940 0.347005    
## temp           -0.07546   0.92731  0.07139 -1.057 0.290488    
## onset_prev_lag -0.36945   0.69111  0.08959 -4.124 3.73e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## airport_dist      0.9126     1.0958    0.7949    1.0477
## conservative      0.9226     1.0838    0.7886    1.0794
## male              0.8700     1.1494    0.7518    1.0068
## age               0.7775     1.2862    0.6503    0.9296
## popdens           0.9125     1.0959    0.7891    1.0553
## manufact          1.0086     0.9914    0.9087    1.1195
## tourism           0.8897     1.1240    0.7912    1.0005
## academics         1.0755     0.9298    0.9174    1.2609
## medinc            1.3456     0.7432    1.1526    1.5709
## healthcare        0.9368     1.0674    0.8177    1.0733
## temp              0.9273     1.0784    0.8062    1.0666
## onset_prev_lag    0.6911     1.4469    0.5798    0.8238
## 
## Concordance= 0.755  (se = 0.015 )
## Likelihood ratio test= 177.8  on 12 df,   p=<2e-16
## Wald test            = 183.9  on 12 df,   p=<2e-16
## Score (logrank) test = 193.2  on 12 df,   p=<2e-16

Predict Slopes

# fixed time windows
lm_slope_prev <- spec_calculate(df = df_ger_prev_scaled, 
               y = "slope_prev", 
               model = "lm", 
               controls = covariates %>% append('slope_prev_lag'),
               all.comb = T)

calc_all_summaries(lm_slope_prev)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04770701 0.047105679 0.9732636 0.4496175 -0.04492057 0.14033459
## sd   0.04866155 0.003961232 0.9445598 0.3249956  0.04463451 0.05352556
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3762331         0.3647529 0.79419554      34.39651 0.0000103227
## sd       0.1110695         0.1112149 0.06709402      12.97273 0.0001987136
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.11260 908.22521 943.66303    235.78389      371.000000
## sd   1.732262   31.40818  60.85522  57.41092     41.98425        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1411133 0.8693848 0.1306152            0.1411133
## sd          0   0.3481809 0.3370202 0.3370202            0.3481809
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.02472864 0.042382347 -0.5866914 0.4324615 -0.10806835 0.05861107
## sd    0.04143443 0.003375866  0.9654986 0.2623563  0.04195917 0.04196624
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean     0.3747063         0.3632077 0.79495961      34.26618 0.001065338
## sd       0.1157000         0.1159113 0.06951753      13.36336 0.018659191
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.39270 908.78539 944.22322    236.36101      371.000000
## sd   1.732262   32.35788  62.76145  59.31843     43.73461        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.09521484 0.3103027 0.6896973                    0
## sd          0  0.29354730 0.4626740 0.4626740                    0
##      significant_negative
## mean           0.09521484
## sd             0.29354730
## 
## $pers_e
##        estimate   std.error statistic   p.value   conf.low  conf.high
## mean 0.09475842 0.042916475  2.145207 0.1152644 0.01036841 0.17914843
## sd   0.05415151 0.002851794  1.072831 0.1384433 0.04912212 0.05928462
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3835310         0.3721721 0.78983225      35.33176 6.972021e-08
## sd       0.1038293         0.1038635 0.06319739      12.45035 1.032508e-06
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -443.15352 904.30703 939.74486    233.02529      371.000000
## sd   1.732262   29.85157  57.72422  54.26548     39.24747        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.4655762        1        0            0.4655762
## sd          0   0.4988745        0        0            0.4988745
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04444294 0.042531589 1.0138143 0.4314829 -0.03919023 0.12807612
## sd   0.03737716 0.003433074 0.8010556 0.3214330  0.03396617 0.04161167
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean     0.3757420         0.3642556 0.7944327      34.36378 0.0001083107
## sd       0.1128081         0.1129825 0.0679861      13.13045 0.0019291217
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -445.19538 908.3908 943.82859    235.96953      371.000000
## sd   1.732262   31.75005  61.5464  58.11313     42.64145        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379   0.1437988 0.9528809 0.04711914            0.1437988
## sd          0   0.3509285 0.2119195 0.21191954            0.3509285
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error statistic   p.value    conf.low  conf.high
## mean -0.06378630 0.042452946 -1.458345 0.3027720 -0.14726483 0.01969223
## sd    0.04686406 0.003144902  0.979003 0.2759218  0.05132561 0.04283259
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3785076         0.3670643 0.79281621      34.71477 2.890617e-06
## sd       0.1093210         0.1094544 0.06617581      12.90910 5.052404e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -444.48386 906.96772 942.40555    234.92414      371.000000
## sd   1.732262   31.04716  60.14311  56.72546     41.32336        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.3442383        0        1                    0
## sd          0   0.4751772        0        0                    0
##      significant_negative
## mean            0.3442383
## sd              0.4751772
plot_all_curves(lm_slope_prev, 'lm_slope_prev')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_slope_prev_ctrl <- lm(slope_prev ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + slope_prev_lag, 
                             data = df_ger_prev_scaled)

summary(lm_slope_prev_ctrl)
## 
## Call:
## lm(formula = slope_prev ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + slope_prev_lag, data = df_ger_prev_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5022 -0.4555 -0.1231  0.3315  3.8222 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -0.007574   0.036347  -0.208 0.835056    
## airport_dist   -0.010112   0.046450  -0.218 0.827784    
## conservative    0.007307   0.056396   0.130 0.896981    
## male           -0.085540   0.053273  -1.606 0.109207    
## age            -0.296574   0.069586  -4.262 2.58e-05 ***
## popdens        -0.116212   0.058141  -1.999 0.046368 *  
## manufact        0.061335   0.042127   1.456 0.146268    
## tourism        -0.126881   0.041916  -3.027 0.002645 ** 
## academics       0.025668   0.056486   0.454 0.649799    
## medinc          0.239428   0.052239   4.583 6.29e-06 ***
## healthcare     -0.018543   0.047757  -0.388 0.698034    
## temp           -0.174761   0.049146  -3.556 0.000426 ***
## slope_prev_lag  0.432857   0.060688   7.133 5.31e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7073 on 366 degrees of freedom
## Multiple R-squared:  0.5156, Adjusted R-squared:  0.4997 
## F-statistic: 32.47 on 12 and 366 DF,  p-value: < 2.2e-16
# variable time windows
lm_slope_prev_var <- spec_calculate(df = df_ger_prev_scaled, 
               y = "slope_prev_var", 
               model = "lm", 
               controls = covariates %>% append('slope_prev_var_lag'),
               all.comb = T)

calc_all_summaries(lm_slope_prev_var)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low conf.high
## mean 0.03138877 0.048027782 0.6168286 0.5240564 -0.06305202 0.1258296
## sd   0.04551793 0.003904359 0.8881225 0.3037374  0.04165726 0.0502624
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3513107         0.3393623 0.81020283      30.71510 0.0001379822
## sd       0.1089591         0.1089781 0.06488532      11.75029 0.0027232807
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.79671 923.59341 959.03124    245.20455      371.000000
## sd   1.732262   29.99134  57.97282  54.44501     41.18655        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.07739258 0.7424316 0.2575684           0.07739258
## sd          0  0.26724596 0.4373484 0.4373484           0.26724596
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.005412083 0.043251658 -0.1122928 0.4938687 -0.09046118 0.07963702
## sd    0.041568173 0.003203958  0.9510700 0.2765537  0.04319003 0.04086338
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean     0.3505259         0.3385690 0.81057266      30.65762 0.001475911
## sd       0.1118288         0.1118893 0.06636315      11.99597 0.025347073
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -452.92184 923.8437 959.28150    245.50122      371.000000
## sd   1.732262   30.55817  59.1121  55.58766     42.27129        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.04272461 0.4780273 0.5219727                    0
## sd          0  0.20226024 0.4995780 0.4995780                    0
##      significant_negative
## mean           0.04272461
## sd             0.20226024
## 
## $pers_e
##        estimate   std.error statistic   p.value    conf.low conf.high
## mean 0.06672988 0.043867749  1.464399 0.3014781 -0.01953070 0.1529904
## sd   0.05252357 0.002821869  1.068661 0.2892869  0.04787854 0.0573294
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3555855         0.3437053 0.80774108      31.20591 2.225108e-06
## sd       0.1036970         0.1036339 0.06204795      11.33341 3.180183e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -451.73375 921.46749 956.90532    243.58870      371.000000
## sd   1.732262   28.85099  55.67972  52.14376     39.19745        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379   0.2744141 0.9606934 0.03930664            0.2744141
## sd          0   0.4462730 0.1943472 0.19434724            0.4462730
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05484896 0.043310646 1.2488281 0.3161042 -0.03031614 0.14001406
## sd   0.03452171 0.003320092 0.7388706 0.2662508  0.03219300 0.03784622
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean     0.3528556         0.3409354 0.8092101      30.95031 0.0001543231
## sd       0.1093621         0.1094004 0.0651510      11.88643 0.0025729140
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000  -452.3210 922.64193 958.07975    244.62057      371.000000
## sd   1.732262    30.1185  58.24095  54.74094     41.33889        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1984863 0.9831543 0.0168457            0.1984863
## sd          0   0.3989090 0.1287089 0.1287089            0.3989090
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.01945947 0.043371756 -0.4091114 0.4680019 -0.10474473 0.06582579
## sd    0.04490624 0.003125042  0.9702821 0.2636699  0.04887248 0.04147441
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3511661         0.3392152 0.81027906      30.70800 7.639954e-05
## sd       0.1092723         0.1093028 0.06506678      11.78065 1.248112e-03
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000  -452.8261 923.65224 959.09007    245.25920      371.000000
## sd   1.732262    30.0687  58.13108  54.60971     41.30493        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.08203125 0.4973145 0.5026855                    0
## sd          0  0.27444583 0.5000538 0.5000538                    0
##      significant_negative
## mean           0.08203125
## sd             0.27444583
plot_all_curves(lm_slope_prev_var, 'lm_slope_prev_var')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_slope_prev_var_ctrl <- lm(slope_prev_var ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + slope_prev_var_lag, 
                             data = df_ger_prev_scaled)

summary(lm_slope_prev_var_ctrl)
## 
## Call:
## lm(formula = slope_prev_var ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + slope_prev_var_lag, data = df_ger_prev_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4605 -0.4684 -0.0894  0.3536  3.8629 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.0107049  0.0370087  -0.289 0.772552    
## airport_dist       -0.0006466  0.0472676  -0.014 0.989094    
## conservative        0.0578272  0.0574076   1.007 0.314452    
## male               -0.1063530  0.0542440  -1.961 0.050679 .  
## age                -0.3348173  0.0704594  -4.752 2.90e-06 ***
## popdens            -0.1144899  0.0591642  -1.935 0.053745 .  
## manufact            0.0844884  0.0428457   1.972 0.049371 *  
## tourism            -0.1367335  0.0426814  -3.204 0.001476 ** 
## academics          -0.0251356  0.0576842  -0.436 0.663279    
## medinc              0.2276969  0.0521603   4.365 1.65e-05 ***
## healthcare         -0.0223008  0.0487212  -0.458 0.647423    
## temp               -0.1940683  0.0505286  -3.841 0.000145 ***
## slope_prev_var_lag  0.4329906  0.0602916   7.182 3.88e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7199 on 366 degrees of freedom
## Multiple R-squared:  0.4982, Adjusted R-squared:  0.4817 
## F-statistic: 30.28 on 12 and 366 DF,  p-value: < 2.2e-16

Predict Cases

lm_cases_0331 <- spec_calculate(df = df_ger_death_scaled, 
               y = "case_rate_0331", 
               model = "lm", 
               controls = covariates %>% append('case_rate_0331_lag'),
               all.comb = T)

calc_all_summaries(lm_cases_0331)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.03158699 0.051022780 0.5994671 0.5627758 -0.06874312 0.13191710
## sd   0.04075752 0.003039436 0.7655479 0.2958935  0.03824516 0.04394422
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27125522         0.2577158 0.86029794     20.714339 0.0002426819
## sd      0.08294606         0.0823706 0.04660655      7.515557 0.0038974652
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.16370 970.32739 1005.76522    275.46553      371.000000
## sd   1.732262   20.54957  39.18622   36.10197     31.35361        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.06103516 0.7888184 0.2111816           0.06103516
## sd          0  0.23942402 0.4081968 0.4081968           0.23942402
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.004091242 0.045901847 -0.08596589 0.5796038 -0.09435164 0.08616916
## sd    0.034700872 0.002269187  0.74910889 0.2613804  0.03536478 0.03460409
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean    0.27028449        0.25673269 0.8607943     20.626725 0.002355802
## sd      0.08560184        0.08506666 0.0479699      7.724115 0.031524022
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -476.35290 970.70581 1006.14364     275.8325      371.000000
## sd   1.732262   21.07044  40.23365   37.14286      32.3575        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.4645996 0.5354004                    0
## sd          0           0 0.4988061 0.4988061                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate   std.error statistic   p.value     conf.low  conf.high
## mean 0.08714459 0.046581747 1.8422270 0.1520497 -0.004452756 0.17874193
## sd   0.04647008 0.001763032 0.9112771 0.1493776  0.043558178 0.04945351
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27812366        0.26469914 0.85640497     21.418135 2.574047e-06
## sd      0.07660416        0.07592723 0.04326494      7.115188 3.140562e-05
##        fit_df fit_logLik  fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -474.51363 967.0273 1002.46509    272.86926      371.000000
## sd   1.732262   19.24243  36.5560   33.48926     28.95637        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.3618164        1        0            0.3618164
## sd          0   0.4805847        0        0            0.4805847
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.07275316 0.045946320 1.5732114 0.1733925 -0.01759469 0.16310101
## sd   0.02970329 0.002372876 0.5983738 0.1394271  0.02789450 0.03209363
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27494180        0.26146891 0.85812315     21.120691 0.0001372159
## sd      0.08246331        0.08190365 0.04643572      7.589597 0.0017659528
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.20596 968.41191 1003.8497    274.07200      371.000000
## sd   1.732262   20.51326  39.13057   36.0854     31.17113        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2634277        1        0            0.2634277
## sd          0   0.4405462        0        0            0.4405462
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error statistic   p.value    conf.low  conf.high
## mean -0.06299523 0.045987096 -1.352350 0.2837433 -0.15342327 0.02743280
## sd    0.03825584 0.002107881  0.770264 0.2261862  0.04073718 0.03608117
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.27430368        0.26081464 0.85855349      21.04466 2.747966e-05
## sd      0.08078661        0.08020112 0.04551672       7.43662 3.526301e-04
##        fit_df fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -475.41660 968.83321 1004.27103    274.31321      371.000000
## sd   1.732262   20.13405  38.37129   35.33896     30.53734        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2299805        0        1                    0
## sd          0   0.4208714        0        0                    0
##      significant_negative
## mean            0.2299805
## sd              0.4208714
plot_all_curves(lm_cases_0331, 'lm_cases_0331')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_cases_0331_ctrl <- lm(case_rate_0331 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + case_rate_0331_lag, 
                             data = df_ger_death_scaled)

summary(lm_cases_0331_ctrl)
## 
## Call:
## lm(formula = case_rate_0331 ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + case_rate_0331_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4908 -0.4246 -0.1322  0.1935  5.7140 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.005574   0.041192  -0.135 0.892434    
## airport_dist        0.007239   0.052639   0.138 0.890700    
## conservative        0.021561   0.064001   0.337 0.736398    
## male               -0.079280   0.060385  -1.313 0.190032    
## age                -0.287805   0.078147  -3.683 0.000265 ***
## popdens            -0.094136   0.065815  -1.430 0.153484    
## manufact            0.045203   0.047719   0.947 0.344117    
## tourism            -0.094058   0.047518  -1.979 0.048517 *  
## academics          -0.011510   0.063625  -0.181 0.856548    
## medinc              0.184173   0.059260   3.108 0.002032 ** 
## healthcare         -0.036733   0.054269  -0.677 0.498913    
## temp               -0.180424   0.055771  -3.235 0.001327 ** 
## case_rate_0331_lag  0.404243   0.071769   5.633 3.54e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8017 on 366 degrees of freedom
## Multiple R-squared:  0.3777, Adjusted R-squared:  0.3573 
## F-statistic: 18.51 on 12 and 366 DF,  p-value: < 2.2e-16
lm_cases_0930 <- spec_calculate(df = df_ger_death_scaled, 
               y = "case_rate_0930", 
               model = "lm", 
               controls = covariates %>% append('case_rate_0930_lag'),
               all.comb = T)

calc_all_summaries(lm_cases_0930)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.04703511 0.048759556 0.9328083 0.4213025 -0.04884462 0.14291485
## sd   0.04623303 0.003708065 0.8870397 0.3049813  0.04263197 0.05063419
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3311692         0.3188131 0.82304322      28.08371 0.0000322385
## sd       0.1042954         0.1042904 0.06154394      11.21308 0.0004799432
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -458.89825 935.79650 971.23433    252.81805      371.000000
## sd   1.732262   28.15626  54.43081  51.21137     39.42368        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1430664 0.8413086 0.1586914            0.1430664
## sd          0   0.3501833 0.3654327 0.3654327            0.3501833
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.002610167 0.044018978 -0.02673773 0.5244260 -0.08916812 0.08394778
## sd    0.040295483 0.003040096  0.89565674 0.2664612  0.04373607 0.03749754
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3292839         0.3169004 0.82408293      27.87911 0.0006614882
## sd       0.1074214         0.1074516 0.06314984      11.42535 0.0139028188
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.32868 936.65736 972.09518    253.53069      371.000000
## sd   1.732262   28.77586  55.66624  52.42724     40.60527        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.0324707 0.5195312 0.4804688                    0
## sd          0   0.1772682 0.4996794 0.4996794                    0
##      significant_negative
## mean            0.0324707
## sd              0.1772682
## 
## $pers_e
##        estimate   std.error statistic    p.value   conf.low  conf.high
## mean 0.11873043 0.044310503 2.6373935 0.03352240 0.03159922 0.20586163
## sd   0.04761501 0.002634405 0.9060562 0.04271978 0.04339097 0.05201209
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.34275305        0.33059563 0.81616073      29.48154 1.913204e-08
## sd      0.09612161        0.09600844 0.05733096      10.80161 3.024722e-07
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -455.83276 929.6655 965.10334    248.43935      371.000000
## sd   1.732262   26.55416  51.2221  48.02528     36.33397        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.7509766        1        0            0.7509766
## sd          0   0.4325002        0        0            0.4325002
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean 0.001807263 0.044142706 0.02930601 0.6298433 -0.08499398 0.08860851
## sd   0.029309078 0.003273494 0.63816227 0.2379830  0.02826462 0.03165452
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean      0.328552         0.3161559 0.8245256      27.78624 0.0008432176
## sd        0.107648         0.1076773 0.0632661      11.42428 0.0168326235
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.52910 937.0582 972.49602    253.80733      371.000000
## sd   1.732262   28.81963  55.7506  52.50393     40.69093        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.5290527 0.4709473                    0
## sd          0           0 0.4992162 0.4992162                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##          estimate   std.error   statistic   p.value    conf.low  conf.high
## mean -0.001840958 0.044147329 -0.01261312 0.5589962 -0.08865129 0.08496938
## sd    0.037168325 0.003113077  0.79081236 0.2465463  0.04075222 0.03430969
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean     0.3290563         0.3166673 0.82425385      27.83814 0.0001595196
## sd       0.1066349         0.1066540 0.06276337      11.33910 0.0028612546
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -459.41897 936.83794 972.27576    253.61671      371.000000
## sd   1.732262   28.63422  55.37859  52.13274     40.30798        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379  0.01220703 0.5891113 0.4108887                    0
## sd          0  0.10982242 0.4920552 0.4920552                    0
##      significant_negative
## mean           0.01220703
## sd             0.10982242
plot_all_curves(lm_cases_0930, 'lm_cases_0930')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_cases_0930_ctrl <- lm(case_rate_0930 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + case_rate_0930_lag, 
                             data = df_ger_death_scaled)

summary(lm_cases_0930_ctrl)
## 
## Call:
## lm(formula = case_rate_0930 ~ airport_dist + conservative + male + 
##     age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + case_rate_0930_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6305 -0.4216 -0.1215  0.3040  3.5774 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -0.003935   0.038008  -0.104 0.917604    
## airport_dist       -0.045205   0.048579  -0.931 0.352702    
## conservative        0.102075   0.058827   1.735 0.083551 .  
## male               -0.034565   0.055737  -0.620 0.535548    
## age                -0.350952   0.071996  -4.875 1.63e-06 ***
## popdens             0.152294   0.061336   2.483 0.013478 *  
## manufact            0.062070   0.044049   1.409 0.159651    
## tourism            -0.109885   0.043906  -2.503 0.012760 *  
## academics          -0.217018   0.059681  -3.636 0.000316 ***
## medinc              0.190831   0.056091   3.402 0.000742 ***
## healthcare          0.013822   0.049981   0.277 0.782285    
## temp               -0.156873   0.051285  -3.059 0.002385 ** 
## case_rate_0930_lag  0.479824   0.068459   7.009 1.16e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7399 on 366 degrees of freedom
## Multiple R-squared:   0.47,  Adjusted R-squared:  0.4526 
## F-statistic: 27.05 on 12 and 366 DF,  p-value: < 2.2e-16

Predict Deaths

lm_deaths_0331 <- spec_calculate(df = df_ger_death_scaled, 
               y = "death_rate_0331", 
               model = "lm", 
               controls = covariates %>% append('death_rate_0331_lag'),
               all.comb = T)

calc_all_summaries(lm_deaths_0331)
## $pers_o
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05926579 0.056344799  1.048407 0.3391156 -0.05152945 0.17006104
## sd   0.02649895 0.002223663  0.465681 0.2207861  0.02562865 0.02803226
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11230093        0.09564048 0.95084691      6.817236 0.003741489
## sd      0.03249807        0.03024746 0.01580279      1.831876 0.022251556
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.578479 1047.15696 1082.59478    335.55025      371.000000
## sd   1.732262    6.852099   11.52894    9.74724     12.28427        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379  0.02636719 0.9821777 0.01782227           0.02636719
## sd          0  0.16024428 0.1323212 0.13232123           0.16024428
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate    std.error  statistic   p.value    conf.low  conf.high
## mean -0.02176714 0.0507656910 -0.4272184 0.6489583 -0.12159172 0.07805745
## sd    0.02201730 0.0008507798  0.4327864 0.2340294  0.02246288 0.02169193
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11002111        0.09331891 0.95206101      6.651521  0.00963241
## sd      0.03320389        0.03098238 0.01615476      1.883209  0.05665206
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -515.059950 1048.11990 1083.557726    336.41202      371.000000
## sd   1.732262    6.973763   11.79201    9.998775     12.55107        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1816406 0.8183594                    0
## sd          0           0 0.3855951 0.3855951                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate    std.error statistic   p.value    conf.low  conf.high
## mean 0.06219816 0.0515241068 1.2102704 0.2768198 -0.03911776 0.16351408
## sd   0.02481655 0.0007406798 0.4899089 0.1810235  0.02544552 0.02425882
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11323239        0.09658630 0.95035178      6.898352 0.002436617
## sd      0.03213563        0.02994777 0.01566049      1.848994 0.014602120
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.381924 1046.76385 1082.201674    335.19815      371.000000
## sd   1.732262    6.787753   11.45944    9.857845     12.14727        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379  0.09594727        1        0           0.09594727
## sd          0  0.29455487        0        0           0.29455487
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.07596961 0.050837168 1.4956105 0.1615796 -0.02399553 0.17593475
## sd   0.01854741 0.001137139 0.3694743 0.1092908  0.01880931 0.01855325
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic fit_p.value
## mean    0.11477591        0.09816324 0.9495095      6.996946 0.003652699
## sd      0.03348276        0.03131897 0.0163866      1.927123 0.022385471
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.040955 1046.08191 1081.51974    334.61471      371.000000
## sd   1.732262    7.080364   12.03403   10.26139     12.65648        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.1005859        1        0            0.1005859
## sd          0   0.3008164        0        0            0.3008164
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.05389933 0.050903492 -1.0622345 0.3222672 -0.15399488 0.04619623
## sd    0.01914286 0.001035423  0.3829244 0.1659477  0.01845582 0.02001427
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11220076        0.09553655 0.95090074      6.820286 0.003463388
## sd      0.03252479        0.03033547 0.01585213      1.860572 0.020486709
##        fit_df  fit_logLik    fit_AIC     fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -514.599597 1047.19919 1082.637020    335.58811      371.000000
## sd   1.732262    6.860555   11.59847    9.943315     12.29437        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379 0.005859375        0        1                    0
## sd          0 0.076331286        0        0                    0
##      significant_negative
## mean          0.005859375
## sd            0.076331286
plot_all_curves(lm_deaths_0331, 'lm_deaths_0331')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_deaths_0331_ctrl <- lm(death_rate_0331 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + death_rate_0331_lag, 
                             data = df_ger_death_scaled)

summary(lm_deaths_0331_ctrl)
## 
## Call:
## lm(formula = death_rate_0331 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + death_rate_0331_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8867 -0.4742 -0.2131  0.1644  5.1407 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.003276   0.047481  -0.069 0.945030    
## airport_dist         0.063879   0.060650   1.053 0.292924    
## conservative         0.033146   0.073530   0.451 0.652415    
## male                -0.170458   0.069611  -2.449 0.014805 *  
## age                 -0.282595   0.087552  -3.228 0.001360 ** 
## popdens             -0.056660   0.075956  -0.746 0.456170    
## manufact             0.190338   0.055019   3.460 0.000605 ***
## tourism             -0.071972   0.054777  -1.314 0.189699    
## academics           -0.081969   0.072606  -1.129 0.259654    
## medinc               0.061189   0.065733   0.931 0.352534    
## healthcare          -0.095838   0.062240  -1.540 0.124469    
## temp                -0.120116   0.064268  -1.869 0.062423 .  
## death_rate_0331_lag  0.231548   0.087866   2.635 0.008766 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.924 on 366 degrees of freedom
## Multiple R-squared:  0.1733, Adjusted R-squared:  0.1462 
## F-statistic: 6.392 on 12 and 366 DF,  p-value: 2.503e-10
lm_deaths_0930 <- spec_calculate(df = df_ger_death_scaled, 
               y = "death_rate_0930", 
               model = "lm", 
               controls = covariates %>% append('death_rate_0930_lag'),
               all.comb = T)

calc_all_summaries(lm_deaths_0930)
## $pers_o
##        estimate   std.error statistic  p.value    conf.low  conf.high
## mean 0.09176990 0.055963743 1.6313826 0.145225 -0.01827605 0.20181584
## sd   0.02910357 0.002231346 0.4896016 0.137522  0.02703758 0.03164673
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.12375372         0.1073184 0.94465296      7.626281 0.002057778
## sd      0.03576864         0.0336473 0.01767784      2.103161 0.013217159
##        fit_df  fit_logLik    fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.087716 1042.17543 1077.6133    331.22109      371.000000
## sd   1.732262    7.626127   13.09662   11.0166     13.52055        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.2592773        1        0            0.2592773
## sd          0   0.4382916        0        0            0.4382916
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate    std.error  statistic   p.value    conf.low  conf.high
## mean -0.02373502 0.0505987501 -0.4672378 0.6381103 -0.12323134 0.07576129
## sd    0.02148019 0.0009822833  0.4226164 0.2372276  0.02203167 0.02109180
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11782829        0.10128121 0.94783532      7.197791 0.009895088
## sd      0.03662981        0.03456711 0.01808525      2.168291 0.058343172
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.359819 1044.71964 1080.15746    333.46091      371.000000
## sd   1.732262    7.745196   13.38859   11.39777     13.84607        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1352539 0.8647461                    0
## sd          0           0 0.3420363 0.3420363                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_e
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.05256399 0.051225663 1.0284875 0.3487208 -0.04816508 0.15329305
## sd   0.02270246 0.000942625 0.4474427 0.1962162  0.02322695 0.02232001
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11989229        0.10338102 0.94673451      7.358761 0.003794686
## sd      0.03579107        0.03374005 0.01768661      2.135709 0.021709173
##        fit_df  fit_logLik   fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -512.922157 1043.8443 1079.28214    332.68071      371.000000
## sd   1.732262    7.596441   13.1105   11.22513     13.52902        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379  0.04541016        1        0           0.04541016
## sd          0  0.20822742        0        0           0.20822742
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##        estimate  std.error statistic   p.value    conf.low  conf.high
## mean 0.01004584 0.05075540 0.1959914 0.7221615 -0.08975852 0.10985019
## sd   0.02082310 0.00116779 0.4083976 0.1840059  0.02043765 0.02144881
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11736593         0.1008111 0.94808125      7.161593  0.01013347
## sd      0.03683933         0.0347719 0.01819243      2.173253  0.06061994
##        fit_df  fit_logLik    fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.457345 1044.91469 1080.3525    333.63568      371.000000
## sd   1.732262    7.789187   13.46958   11.4449     13.92527        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.6926270 0.3073730                    0
## sd          0           0 0.4614616 0.4614616                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_n
##        estimate   std.error statistic   p.value    conf.low  conf.high
## mean 0.03733565 0.050718234 0.7335759 0.4914366 -0.06239562 0.13706691
## sd   0.02251156 0.001073727 0.4357044 0.2497344  0.02195929 0.02324317
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic fit_p.value
## mean    0.11858409        0.10205446 0.94741952      7.243279  0.01153319
## sd      0.03739883        0.03533075 0.01848994      2.202537  0.06886308
##        fit_df  fit_logLik    fit_AIC    fit_BIC fit_deviance fit_df.residual
## mean 7.000000 -513.190383 1044.38077 1079.81859    333.17521      371.000000
## sd   1.732262    7.912561   13.70082   11.59235     14.13676        1.732262
##      fit_nobs significant  positive   negative significant_positive
## mean      379           0 0.9448242 0.05517578                    0
## sd          0           0 0.2283509 0.22835092                    0
##      significant_negative
## mean                    0
## sd                      0
plot_all_curves(lm_deaths_0930, 'lm_deaths_0930')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_deaths_0930_ctrl <- lm(death_rate_0930 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + death_rate_0930_lag, 
                             data = df_ger_death_scaled)

summary(lm_deaths_0930_ctrl)
## 
## Call:
## lm(formula = death_rate_0930 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + death_rate_0930_lag, data = df_ger_death_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9553 -0.5092 -0.2248  0.2556  6.5021 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)         -0.004634   0.047233  -0.098  0.92190   
## airport_dist         0.026542   0.060328   0.440  0.66023   
## conservative         0.026796   0.073079   0.367  0.71407   
## male                -0.087034   0.069386  -1.254  0.21052   
## age                 -0.211738   0.086909  -2.436  0.01531 * 
## popdens              0.029944   0.075712   0.396  0.69270   
## manufact             0.135993   0.054799   2.482  0.01353 * 
## tourism             -0.093115   0.054578  -1.706  0.08884 . 
## academics           -0.162248   0.072481  -2.238  0.02579 * 
## medinc               0.093815   0.065408   1.434  0.15234   
## healthcare          -0.026707   0.061800  -0.432  0.66588   
## temp                -0.184670   0.064375  -2.869  0.00436 **
## death_rate_0930_lag  0.284429   0.087564   3.248  0.00127 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9191 on 366 degrees of freedom
## Multiple R-squared:  0.182,  Adjusted R-squared:  0.1552 
## F-statistic: 6.788 on 12 and 366 DF,  p-value: 4.504e-11

Predict Socdist Onset

cox_onset_socdist <- spec_calculate(df = df_ger_socdist_scaled, 
               y = "Surv(cpt_day_socdist_2, event)", 
               model = "cox_model", 
               controls = covariates %>% 
                 append(c('cpt_day_socdist_2_lag', 'onset_prev', 'slope_prev')),
               all.comb = T)

cox_onset_socdist <- cox_onset_socdist %>% map(filter_socdist)

cox_onset_socdist_hr <- cox_onset_socdist %>% map(coef_to_hr)

calc_all_summaries(cox_onset_socdist_hr)
## $pers_o
##        estimate   std.error  statistic   p.value   conf.low  conf.high fit_n
## mean 0.92412115 0.060146732 -1.3205350 0.2201715 0.82135075 1.03977081   379
## sd   0.02156262 0.002248289  0.4019257 0.1387493 0.01899451 0.02532624     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.32341     0.005706523         36.55378
## sd            0          13.75696     0.024395795         12.35949
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.005739844           36.14392      0.005735019                   NA
## sd      0.023985584           12.05654      0.023777053                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09557573         0.9999496      0.61749673
## sd                   NA    0.03296967         0.0000000      0.03336312
##      fit_std.error.concordance   fit_logLik    fit_AIC  fit_BIC fit_nobs
## mean               0.028868295 -1856.052442 3730.10488 3765.543      379
## sd                 0.001210389     6.878482   12.02107   11.250        0
##      significant positive negative significant_positive significant_negative
## mean  0.06811523        0        1                    0           0.06811523
## sd    0.25197430        0        0                    0           0.25197430
## 
## $pers_c
##        estimate   std.error  statistic   p.value  conf.low  conf.high fit_n
## mean 0.92379825 0.053425612 -1.4965299 0.1765796 0.8319259 1.02582137   379
## sd   0.02216087 0.001129165  0.4722299 0.1393372 0.0186606 0.02630592     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.89125     0.002080154         37.09716
## sd            0          12.54155     0.007141225         11.33829
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.002219141           36.66871      0.002244855                   NA
## sd      0.007425369           11.03882      0.007438921                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09703073         0.9999496      0.62825725
## sd                   NA    0.02998817         0.0000000      0.02806221
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean                0.02888723 -1855.768524 3729.53705 3764.97487      379
## sd                  0.00105630     6.270773   10.78844   10.29674        0
##      significant positive negative significant_positive significant_negative
## mean   0.1914062        0        1                    0            0.1914062
## sd     0.3934561        0        0                    0            0.3934561
## 
## $pers_e
##        estimate    std.error  statistic    p.value   conf.low  conf.high fit_n
## mean 0.88750828 0.0562011720 -2.1242714 0.04010842 0.79494811 0.99084875   379
## sd   0.01407696 0.0008904741  0.2702866 0.02493297 0.01325974 0.01506055     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          41.01862     0.002206119         39.25972
## sd            0          13.53017     0.009699107         12.25642
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean     0.00234413           38.82729      0.002322961                   NA
## sd       0.01021985           11.92044      0.010047201                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA     0.1020040         0.9999496      0.62851368
## sd                   NA     0.0321979         0.0000000      0.02755481
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028252338 -1854.704839 3727.40968 3762.84750      379
## sd                 0.001407939     6.765084   11.72497   10.85692        0
##      significant positive negative significant_positive significant_negative
## mean   0.6994629        0        1                    0            0.6994629
## sd     0.4585476        0        0                    0            0.4585476
## 
## $pers_a
##        estimate   std.error  statistic   p.value   conf.low conf.high fit_n
## mean 0.94998162 0.055857806 -0.9252410 0.3798245 0.85145243 1.0599184   379
## sd   0.01686576 0.001206494  0.3290921 0.1672532 0.01424408 0.0201586     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          37.40628      0.00487520         35.73613
## sd            0          13.19216      0.01728311         11.95547
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.005074892           35.31908      0.005112128                   NA
## sd      0.017780079           11.64557      0.017763039                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09343289         0.9999496      0.62513121
## sd                   NA    0.03167007         0.0000000      0.03199302
##      fit_std.error.concordance  fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028872705 -1856.51101 3731.02201 3766.45984      379
## sd                 0.001126387     6.59608   11.46145   10.84048        0
##      significant positive negative significant_positive significant_negative
## mean           0        0        1                    0                    0
## sd             0        0        0                    0                    0
## 
## $pers_n
##        estimate   std.error statistic    p.value   conf.low  conf.high fit_n
## mean 1.08916744 0.057407948 1.4897329 0.14528486 0.97327553 1.21886580   379
## sd   0.01220407 0.001200459 0.2138959 0.06005665 0.01235991 0.01223976     0
##      fit_nevent fit_statistic.log fit_p.value.log fit_statistic.sc
## mean        379          38.68176      0.00353503         37.02830
## sd            0          13.43207      0.01343931         12.20186
##      fit_p.value.sc fit_statistic.wald fit_p.value.wald fit_statistic.robust
## mean    0.003639324           36.63856      0.003605216                   NA
## sd      0.013488312           11.86728      0.013282388                   NA
##      fit_p.value.robust fit_r.squared fit_r.squared.max fit_concordance
## mean                 NA    0.09645849         0.9999496      0.62898123
## sd                   NA    0.03214506         0.0000000      0.03045481
##      fit_std.error.concordance   fit_logLik    fit_AIC    fit_BIC fit_nobs
## mean               0.028670041 -1855.873268 3729.74654 3765.18436      379
## sd                 0.001247607     6.716033   11.72373   11.09026        0
##      significant positive negative significant_positive significant_negative
## mean           0        1        0                    0                    0
## sd             0        0        0                    0                    0
plot_all_curves(cox_onset_socdist_hr, 'cox_onset_socdist_hr', hr = T)
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

cox_onset_socdist_ctrl <- coxph(Surv(cpt_day_socdist_2, event) ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev + slope_prev +
                                 cpt_day_socdist_2_lag,
                             data = df_ger_socdist_scaled)

summary(cox_onset_socdist_ctrl)
## Call:
## coxph(formula = Surv(cpt_day_socdist_2, event) ~ airport_dist + 
##     conservative + male + age + popdens + manufact + tourism + 
##     academics + medinc + healthcare + temp + onset_prev + slope_prev + 
##     cpt_day_socdist_2_lag, data = df_ger_socdist_scaled)
## 
##   n= 379, number of events= 379 
## 
##                            coef exp(coef)  se(coef)      z Pr(>|z|)    
## airport_dist           0.125877  1.134142  0.067907  1.854   0.0638 .  
## conservative          -0.331611  0.717767  0.082851 -4.002 6.27e-05 ***
## male                  -0.141864  0.867740  0.081693 -1.737   0.0825 .  
## age                   -0.017427  0.982724  0.103089 -0.169   0.8658    
## popdens                0.108749  1.114883  0.091485  1.189   0.2346    
## manufact              -0.072176  0.930367  0.052325 -1.379   0.1678    
## tourism                0.009962  1.010012  0.061504  0.162   0.8713    
## academics             -0.058010  0.943640  0.079147 -0.733   0.4636    
## medinc                -0.104273  0.900979  0.071349 -1.461   0.1439    
## healthcare            -0.173877  0.840400  0.069290 -2.509   0.0121 *  
## temp                   0.077350  1.080420  0.072919  1.061   0.2888    
## onset_prev            -0.013565  0.986527  0.075247 -0.180   0.8569    
## slope_prev             0.146099  1.157311  0.086200  1.695   0.0901 .  
## cpt_day_socdist_2_lag -0.238041  0.788170  0.114252 -2.083   0.0372 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## airport_dist             1.1341     0.8817    0.9928    1.2956
## conservative             0.7178     1.3932    0.6102    0.8443
## male                     0.8677     1.1524    0.7394    1.0184
## age                      0.9827     1.0176    0.8029    1.2028
## popdens                  1.1149     0.8970    0.9319    1.3338
## manufact                 0.9304     1.0748    0.8397    1.0308
## tourism                  1.0100     0.9901    0.8953    1.1394
## academics                0.9436     1.0597    0.8080    1.1020
## medinc                   0.9010     1.1099    0.7834    1.0362
## healthcare               0.8404     1.1899    0.7337    0.9626
## temp                     1.0804     0.9256    0.9365    1.2464
## onset_prev               0.9865     1.0137    0.8513    1.1433
## slope_prev               1.1573     0.8641    0.9774    1.3703
## cpt_day_socdist_2_lag    0.7882     1.2688    0.6300    0.9860
## 
## Concordance= 0.678  (se = 0.029 )
## Likelihood ratio test= 61.37  on 14 df,   p=7e-08
## Wald test            = 57.4  on 14 df,   p=3e-07
## Score (logrank) test = 58.42  on 14 df,   p=2e-07

Predict Socdist Mean

lm_mean_socdist <- spec_calculate(df = df_ger_socdist_scaled, 
               y = "mean_diff_socdist_2", 
               model = "lm", 
               controls = covariates %>% 
                 append(c('mean_diff_socdist_2_lag', 'onset_prev', 'slope_prev')), 
               all.comb = T)

lm_mean_socdist <- lm_mean_socdist %>% map(filter_socdist)

calc_all_summaries(lm_mean_socdist)
## $pers_o
##        estimate   std.error statistic    p.value   conf.low  conf.high
## mean 0.15330885 0.044355910  3.425476 0.02297009 0.06608679 0.24053090
## sd   0.07755911 0.001414087  1.656772 0.03632131 0.07585609 0.07932302
##      fit_r.squared fit_adj.r.squared fit_sigma fit_statistic  fit_p.value
## mean    0.45206387        0.43888864 0.7476768     35.209876 7.055131e-20
## sd      0.07030338        0.07031662 0.0457296      7.949236 2.044565e-18
##        fit_df fit_logLik   fit_AIC  fit_BIC fit_deviance fit_df.residual
## mean 9.000000  -421.8148 865.62962 908.9425    207.11986      369.000000
## sd   1.732262    23.1611  44.31669  40.9420     26.57468        1.732262
##      fit_nobs significant positive negative significant_positive
## mean      379   0.8122559        1        0            0.8122559
## sd          0   0.3905554        0        0            0.3905554
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_c
##         estimate  std.error  statistic   p.value    conf.low  conf.high
## mean -0.03851258 0.04069143 -0.9421542 0.4218812 -0.11852873 0.04150356
## sd    0.02790481 0.00308532  0.6560585 0.2885431  0.02975273 0.02730813
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.43015918        0.41651193 0.76145621     32.660794 8.760185e-07
## sd      0.09616314        0.09669053 0.06060868      9.707367 2.157079e-05
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.28432 878.56865 921.88154    215.39983      369.000000
## sd   1.732262   29.50691  57.06175  53.66655     36.34967        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.0769043 0.0546875 0.9453125                    0
## sd          0   0.2664721 0.2273970 0.2273970                    0
##      significant_negative
## mean            0.0769043
## sd              0.2664721
## 
## $pers_e
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.01503002 0.041886341 -0.3896659 0.5811354 -0.09739585 0.06733580
## sd    0.02278655 0.003250974  0.5237195 0.1998854  0.01851793 0.02787933
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.42867767        0.41499876 0.76240500      32.47312 3.773528e-06
## sd      0.09720623        0.09774027 0.06116183       9.72353 8.896818e-05
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.73869 879.4774 922.79028    215.95984      369.000000
## sd   1.732262   29.73164  57.4969  54.06717     36.74395        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.1518555 0.8481445                    0
## sd          0           0 0.3589246 0.3589246                    0
##      significant_negative
## mean                    0
## sd                      0
## 
## $pers_a
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.02332656 0.040841977 -0.5897552 0.4122219 -0.10363874 0.05698563
## sd    0.04379788 0.002890095  1.0885441 0.2944594  0.04279858 0.04549033
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.43026173        0.41662073 0.76130558     32.712822 6.230705e-06
## sd      0.09788511        0.09844834 0.06159559      9.858796 1.500141e-04
##        fit_df fit_logLik  fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.17439 878.3488 921.66168    215.36107      369.000000
## sd   1.732262   29.92968  57.9119  54.51971     37.00057        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379   0.1582031 0.3818359 0.6181641                    0
## sd          0   0.3649759 0.4858960 0.4858960                    0
##      significant_negative
## mean            0.1582031
## sd              0.3649759
## 
## $pers_n
##         estimate   std.error  statistic   p.value    conf.low  conf.high
## mean -0.00604797 0.041254131 -0.1365582 0.7023933 -0.08717062 0.07507468
## sd    0.01968351 0.003240101  0.4720721 0.2018550  0.02250371 0.01869854
##      fit_r.squared fit_adj.r.squared  fit_sigma fit_statistic  fit_p.value
## mean    0.42841045        0.41472473 0.76258917     32.435635 6.023785e-06
## sd      0.09718258        0.09771424 0.06110589      9.710442 1.448236e-04
##        fit_df fit_logLik   fit_AIC   fit_BIC fit_deviance fit_df.residual
## mean 9.000000 -428.83354 879.66708 922.97998    216.06085      369.000000
## sd   1.732262   29.68545  57.40582  53.97958     36.73502        1.732262
##      fit_nobs significant  positive  negative significant_positive
## mean      379           0 0.4108887 0.5891113                    0
## sd          0           0 0.4920552 0.4920552                    0
##      significant_negative
## mean                    0
## sd                      0
plot_all_curves(lm_mean_socdist, 'lm_mean_socdist')
## $spec_results_o

## 
## $spec_results_c

## 
## $spec_results_e

## 
## $spec_results_a

## 
## $spec_results_n

lm_mean_socdist_ctrl <- lm(mean_diff_socdist_2 ~ airport_dist + 
                               conservative + male + age + popdens + manufact +
                               tourism + academics + medinc + healthcare + 
                               temp + onset_prev + slope_prev +
                                 mean_diff_socdist_2_lag,
                             data = df_ger_socdist_scaled)

summary(lm_mean_socdist_ctrl)
## 
## Call:
## lm(formula = mean_diff_socdist_2 ~ airport_dist + conservative + 
##     male + age + popdens + manufact + tourism + academics + medinc + 
##     healthcare + temp + onset_prev + slope_prev + mean_diff_socdist_2_lag, 
##     data = df_ger_socdist_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.85241 -0.44592  0.01798  0.37505  2.81498 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -0.001481   0.035208  -0.042   0.9665    
## airport_dist            -0.010038   0.045187  -0.222   0.8243    
## conservative            -0.273373   0.055239  -4.949 1.14e-06 ***
## male                    -0.117690   0.051800  -2.272   0.0237 *  
## age                      0.003341   0.067934   0.049   0.9608    
## popdens                  0.328194   0.056804   5.778 1.63e-08 ***
## manufact                 0.085452   0.041025   2.083   0.0380 *  
## tourism                 -0.013999   0.041163  -0.340   0.7340    
## academics                0.364853   0.054490   6.696 8.15e-11 ***
## medinc                  -0.073709   0.050032  -1.473   0.1416    
## healthcare              -0.008339   0.046183  -0.181   0.8568    
## temp                    -0.081000   0.048476  -1.671   0.0956 .  
## onset_prev               0.008089   0.051930   0.156   0.8763    
## slope_prev               0.031435   0.055945   0.562   0.5745    
## mean_diff_socdist_2_lag  0.135695   0.064960   2.089   0.0374 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6853 on 364 degrees of freedom
## Multiple R-squared:  0.5478, Adjusted R-squared:  0.5304 
## F-statistic: 31.49 on 14 and 364 DF,  p-value: < 2.2e-16

Save all results

results_ger <- list(cox_onset_prev_hr, lm_slope_prev, lm_slope_prev_var, 
                   lm_cases_0331, lm_cases_0930, lm_deaths_0331, lm_deaths_0930, 
                   cox_onset_socdist_hr, lm_mean_socdist)


names(results_ger) <- list('cox_onset_prev_hr', 'lm_slope_prev', 'lm_slope_prev_var', 
                   'lm_cases_0331', 'lm_cases_0930', 'lm_deaths_0331', 'lm_deaths_0930', 
                   'cox_onset_socdist_hr', 'lm_mean_socdist')

save(results_ger, file="GER_spec_results.RData")

Descriptives

Distributions

df_ger_desc <- df_ger_socdist %>% 
  dplyr::select(pers_o, pers_c, pers_e, pers_a, pers_n, 
                age, male, conservative, 
                academics, medinc, manufact, 
                airport_dist, tourism, healthcare, popdens) %>% 
  distinct() %>% drop_na()

ger_means <- df_ger_desc %>% summarise_all(mean)
ger_sd <- df_ger_desc %>% summarise_all(sd)

desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame()

desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger  

Explore correlations

a <- df_ger_socdist_unscaled %>% select(kreis, pers_o:pers_n, cpt_day_socdist, mean_diff_socdist)
b <- df_ger_prev_unscaled %>% select(kreis, onset_prev, slope_prev)
c <- df_ger_death_unscaled %>% select(kreis, case_rate_0331:death_rate_0930)

ger_joined <- plyr::join_all(list(a, b, c), by='kreis')

ger_joined %>% select(-kreis) %>% cor(use = 'pairwise.complete')
##                        pers_o      pers_c      pers_e      pers_a        pers_n
## pers_o             1.00000000 -0.10201270  0.27983871  0.15759347 -0.0640347187
## pers_c            -0.10201270  1.00000000  0.21472242  0.32690021 -0.3551086811
## pers_e             0.27983871  0.21472242  1.00000000  0.18288787 -0.3649780097
## pers_a             0.15759347  0.32690021  0.18288787  1.00000000 -0.4280315714
## pers_n            -0.06403472 -0.35510868 -0.36497801 -0.42803157  1.0000000000
## cpt_day_socdist   -0.10416238  0.10437674 -0.04952867 -0.01672102  0.0138442293
## mean_diff_socdist  0.31429946 -0.12682312  0.09840937 -0.08341683  0.0006497013
## onset_prev        -0.16574985  0.05702967 -0.27615283 -0.09470167  0.2352513229
## slope_prev         0.15470785 -0.04204164  0.24306903  0.13157797 -0.1983834025
## case_rate_0331     0.10695243 -0.01752091  0.20572870  0.13949525 -0.1722746724
## death_rate_0331    0.06566400 -0.03556727  0.10805800  0.09283184 -0.0999658180
## case_rate_0930     0.15634007 -0.05283919  0.25572247  0.04675307 -0.1116940920
## death_rate_0930    0.07587651 -0.03740828  0.08958204  0.02599047 -0.0095583322
##                   cpt_day_socdist mean_diff_socdist  onset_prev  slope_prev
## pers_o              -0.1041623824      0.3142994608 -0.16574985  0.15470785
## pers_c               0.1043767434     -0.1268231184  0.05702967 -0.04204164
## pers_e              -0.0495286726      0.0984093720 -0.27615283  0.24306903
## pers_a              -0.0167210202     -0.0834168254 -0.09470167  0.13157797
## pers_n               0.0138442293      0.0006497013  0.23525132 -0.19838340
## cpt_day_socdist      1.0000000000      0.0358262655  0.14392128 -0.11996515
## mean_diff_socdist    0.0358262655      1.0000000000 -0.20519894  0.22948985
## onset_prev           0.1439212762     -0.2051989385  1.00000000 -0.67617865
## slope_prev          -0.1199651501      0.2294898547 -0.67617865  1.00000000
## case_rate_0331      -0.0689721829      0.1733238658 -0.66024501  0.92767374
## death_rate_0331      0.0003684063      0.1435286041 -0.40795429  0.69587384
## case_rate_0930      -0.1260015001      0.2185131719 -0.58477976  0.82060844
## death_rate_0930     -0.0078242968      0.1456354879 -0.32697005  0.64249969
##                   case_rate_0331 death_rate_0331 case_rate_0930 death_rate_0930
## pers_o                0.10695243    0.0656640020     0.15634007     0.075876508
## pers_c               -0.01752091   -0.0355672682    -0.05283919    -0.037408282
## pers_e                0.20572870    0.1080579960     0.25572247     0.089582044
## pers_a                0.13949525    0.0928318364     0.04675307     0.025990472
## pers_n               -0.17227467   -0.0999658180    -0.11169409    -0.009558332
## cpt_day_socdist      -0.06897218    0.0003684063    -0.12600150    -0.007824297
## mean_diff_socdist     0.17332387    0.1435286041     0.21851317     0.145635488
## onset_prev           -0.66024501   -0.4079542869    -0.58477976    -0.326970047
## slope_prev            0.92767374    0.6958738421     0.82060844     0.642499692
## case_rate_0331        1.00000000    0.7967609328     0.82553862     0.687015730
## death_rate_0331       0.79676093    1.0000000000     0.68600498     0.850162063
## case_rate_0930        0.82553862    0.6860049833     1.00000000     0.745023167
## death_rate_0930       0.68701573    0.8501620632     0.74502317     1.000000000
ger_means <- ger_joined %>% select(-kreis) %>% summarise_all(mean)
ger_sd <- ger_joined %>% select(-kreis) %>% summarise_all(sd)

desc_ger <- rbind(ger_means, ger_sd) %>% t() %>% round(3) %>% as.data.frame() %>%
  rename(mean=V1, sd=V2)

desc_ger %>% rownames_to_column() %>% write_csv('ger_descriptives.csv')
desc_ger  

Visualizations

df_ger_prev <- read_csv('GER_prevalence.csv')

df_ger_prev <- df_ger_prev %>%
  mutate(date = as.Date(date, "%d%b%Y"),
         kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day) %>%
  filter(date < '2020-04-01')

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence))
names(df_dates) <- c('date', 'time')

tto <- df_ger_prev_unscaled %>%
  filter(event==1) %>%
  merge(df_dates, by.x = 'onset_prev', by.y = 'time') %>%
  ggplot(aes(x=date+1)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - COVID-19 Onset') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Date')
  
soa <- df_ger_slope_prev %>% 
  ggplot(aes(x=slope_prev)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - COVID-19 Growth Rates') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Standardized growth rates')

figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)

figure

tto <- df_ger_death %>% 
  ggplot(aes(x=case_rate_0930)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Cumulative Case Rates') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Cases Per 1000 Inhabitants')
  
soa <- df_ger_death %>%
  ggplot(aes(x=death_rate_0930)) +
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Cumulative Death Rates') +
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Deaths Per 1000 Inhabitants')


figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)
figure

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = data.frame(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')


tto <- df_ger_cpt_socdist %>% 
  merge(df_dates, by.x = 'cpt_day_socdist_2', by.y = 'time') %>%
  #filter(cpt_day_socdist_2 < 30) %>% 
  ggplot(aes(x=date+1)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Time to Adoption') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Date') +
  xlim(c(as.Date('2020-03-01'), as.Date('2020-03-31')))
  
soa <- df_ger_cpt_socdist %>% 
  ggplot(aes(x=mean_diff_socdist_2)) + 
  geom_histogram(bins=15) +
  theme_light() +
  ggtitle('Distribution - Strength of Adjustment') + 
  theme(plot.title = element_text(size=11)) +
  ylab('Absolute frequency') +
  xlab('Standardized mean difference')


figure <- ggarrange(tto, soa,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)
figure

>>>>>>> f9bb205265c8cd2a54b17b2b2634fb89e3d3c38f